module ncepnems_io !$$$ module documentation block ! . . . . ! module: ncepnems_io ! prgmmr: Huang org: np23 date: 2010-02-22 ! ! abstract: This module contains routines which handle input/output ! operations for NCEP NEMS global atmospheric and surface files. ! ! program history log: ! 2010-02-22 Huang Initial version. Based on ncepgfs_io ! 2010-10-18 Huang Remove subroutine reorder_gfsgrib for no longer been called in GSI ! For Now, subroutine sfc_interpolate is kept in ncepgfs_io.f90. ! When sigio and gfsio are both retired, i.e., remove ncepgfs_io.f90. ! move this routines back to this module ! 2011-03-03 Huang Changes has been made to adopt to high resolution GSI run (T382 & T574) ! both for CPU and memory issues. ! Future development of nemsio need to consider a mapping routine be ! inserted between read-in forecast field and GSI first guess field, ! as well as GSI analysis field and write-out data field for forecast ! model. Due to computation resource, GSI may not be able to run at ! the same resolution as that of forecast model, e.g., run GSI at T382 ! w/ T574 forecast model output. ! 2011-10-25 Huang (1) Add unified error message routine to make the code cleaner ! (2) To reduce the memory allocation as low as possible, remove all ! reference to sfc_head and re-used the same local arrays. ! Remove unneeded nemsio_data & gfsdata. ! (3) Add parallel IO code to read_atm_ ! 2011-11-01 Huang (1) add integer nst_gsi to control the mode of NSST ! (2) add read_nemsnst to read ncep nst file ! (3) add subroutine write_nemssfc_nst to save sfc and nst files ! 2016-01-01 Li (1) Move tran_gfssfc from ncepgfs_io.f90 to here ! (2) Modify write_sfc_nst_ to follows the update done in sfcio ! (3) Modify read_sfc_ to follows the update done in sfcio for more effective I/O ! 2016-04-20 Li Modify to handle the updated nemsio sig file (P, DP & DPDT removed) ! 2016-08-18 li - tic591: add read_sfc_anl & read_nemssfc_anl to read nemsio sfc file (isli only) with analysis resolution ! change/modify sfc_interpolate to be intrp22 to handle more general interpolation (2d to 2d) ! 2016-11-18 li - tic615: change nst mask name from slmsk to land ! 2017-08-30 li - tic659: modify read_nems_sfc_ and read_sfc_ to read sfc file in ! nemsio Gaussin grids generated by FV3 WriteComponent ! 2018-05-19 eliu - add I/O component for fv3 hydrometeors ! ! Subroutines Included: ! sub read_nems - driver to read ncep nems atmospheric and surface ! sub read_nems_chem ! sub read_nemsatm - read ncep nems atmospheric file, scatter ! on grid to analysis subdomains ! sub read_nemssfc - read ncep nems surface file, scatter on grid to ! analysis subdomains ! sub read_nemssfc_anl- read ncep EnKF nems surface file, scatter on grid to ! analysis subdomains ! sub write_nems - driver to write ncep nems atmospheric and surface ! analysis files ! sub write_nemsatm - gather on grid, write ncep nems atmospheric analysis file ! sub write_nemssfc - gather/write on grid ncep surface analysis file ! sub read_nemsnst - read ncep nst file, scatter on grid to analysis subdomains ! sub write_nems_sfc_nst - gather/write on grid ncep surface & nst analysis file ! sub intrp22 - interpolate from one grid to another grid (2D) ! sub read_nems_sfcnst - read sfc hist file, including sfc and nst vars, scatter on grid to analysis subdomains ! ! Variable Definitions: ! The difference of time Info between operational GFS IO (gfshead%, sfc_head%), ! analysis time (iadate), and NEMSIO (idate=) ! ! gfshead & sfc_head NEMSIO Header Analysis time (obsmod) ! =================== ============================ ========================== ! %idate(1) Hour idate(1) Year iadate(1) Year ! %idate(2) Month idate(2) Month iadate(2) Month ! %idate(3) Day idate(3) Day iadate(3) Day ! %idate(4) Year idate(4) Hour iadate(4) Hour ! idate(5) Minute iadate(5) Minute ! idate(6) Scaled seconds ! idate(7) Seconds multiplier ! ! The difference of header forecasting hour Info bewteen operational GFS IO ! (gfshead%, sfc_head%) and NEMSIO ! ! gfshead & sfc_head NEMSIO Header ! ========================== ============================ ! %fhour FCST Hour (r_kind) nfhour FCST Hour (i_kind) ! nfminute FCST Mins (i_kind) ! nfsecondn FCST Secs (i_kind) numerator ! nfsecondd FCST Secs (i_kind) denominator ! ! %fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600 ! ! nframe - nframe is the number of grids extend outward from the ! edge of modeling domain. ! ! NEMSIO provides a more flexible read. User can get the ! size of record (1D) to be read from file header. The ! normal record size should be delx*dely, i.e., total model ! grid points. However, some regional models also ouput ! additional data of grids around the modeling domain ! (buffer zone). For this type of output, nframe needs to ! be know to calculate the size of record, i.e., ! array size = (delx+2*nframe) * (dely+2*nframe) ! ! However, nframe should always be zero for global model. ! To simplify the code for reading and writing global model ! files, we will not factor in the nframe for computing ! array size or array index shift (by nframe) between ! input/output array and internal GSI array. The normal ! size of I/O record remains as delx*dely. Add a checking ! routine to assure nframe=zero. ! ! def imp_physics - type of microphysics used in the GFS. 99: Zhao-Carr, 11: GFDL ! def lupp - if T, UPP is used and additional variables are output ! ! attributes: ! language: f90 ! machine: ! ! NOTE: When global meteorology switched to NEMS/GFS, all routines and ! modules of old GFS (sigio) can be deactivated. To keep the code ! clean, all "nems" can be replaced by "gfs" for minimal changes ! of GSI code structure. For dual purpose, two distincit routine ! names are used to accomodiate old and new systems. It is now ! controled by a namelist argument "use_gfs_nemsio" ! ! !$$$ end documentation block use constants, only: zero,one,fv,r60,r3600 use kinds, only: i_kind,r_kind implicit none private public init_nems public read_nems public read_nems_chem public read_nemsatm public read_nemssfc public read_nemssfc_anl public berror_read_hsst public write_fv3atm_nems public write_nemsatm public write_nemssfc public read_nemsnst public write_nems_sfc_nst public intrp22 public tran_gfssfc public error_msg public imp_physics public lupp interface init_nems module procedure init_ end interface interface read_nems module procedure read_ end interface interface read_nems_chem module procedure read_chem_ end interface interface read_nemsatm module procedure read_atm_ end interface interface read_nemssfc module procedure read_nemssfc_ end interface interface read_nemssfc_anl module procedure read_nemssfc_anl_ end interface interface berror_read_hsst module procedure berror_read_hsst_ end interface interface read_nemsnst module procedure read_nemsnst_ end interface interface write_nemsatm module procedure write_atm_ end interface interface write_fv3atm_nems module procedure write_fv3atm_ end interface interface write_nemssfc module procedure write_sfc_ end interface interface write_nems_sfc_nst module procedure write_sfc_nst_ end interface interface error_msg module procedure error_msg_ end interface character(len=*),parameter::myname='ncepnems_io' integer(i_kind),save :: imp_physics logical,save :: lupp contains subroutine init_ !$$$ subprogram documentation block ! . . . ! subprogram: read_nems ! ! prgrmmr: Todling ! ! abstract: ! ! program history log: ! 2019-07-11 Todling - create to initialize vars that should not be in CV file ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block imp_physics=99 lupp = .false. end subroutine init_ subroutine read_ !$$$ subprogram documentation block ! . . . ! subprogram: read_nems ! ! prgrmmr: Ho-Chun Huang ! ! abstract: ! ! program history log: ! 2010-03-31 Huang - create routine based on read_gfs ! 2010-10-19 Huang - remove spectral part for gridded NEMS/GFS ! 2011-05-01 todling - cwmr no longer in guess-grids; use metguess bundle now ! 2013-10-19 todling - metguess now holds background ! 2016-03-30 todling - update interface to general read (pass bundle) ! 2016-06-23 Li - Add cloud partitioning, which was missed (based on GFS ! ticket #239, comment 18) ! 2018-05-19 eliu - add components to read in hydrometeor related ! variables ! 2019-07-10 Zhu - Add convective clouds ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: i_kind,r_kind use constants, only: qcmin use gridmod, only: sp_a,grd_a,lat2,lon2,nsig use guess_grids, only: ifilesig,nfldsig,ntguessig use gsi_metguess_mod, only: gsi_metguess_bundle,gsi_metguess_get use guess_grids, only: ifilesig,nfldsig use gsi_metguess_mod, only: gsi_metguess_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use gsi_bundlemod, only: gsi_bundlecreate use gsi_bundlemod, only: gsi_grid use gsi_bundlemod, only: gsi_gridcreate use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundledestroy use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info use mpimod, only: npe,mype use cloud_efr_mod, only: cloud_calc_gfs,set_cloud_lower_bound use gridmod, only: fv3_full_hydro implicit none character(len=*),parameter::myname_=myname//'*read_' character(24) filename integer(i_kind):: it, istatus, inner_vars, num_fields integer(i_kind):: i,j,k real(r_kind),pointer,dimension(:,: ):: ges_ps_it =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_z_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_u_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_v_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_div_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_vor_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_tv_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_q_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_oz_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_cwmr_it=>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_ql_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qi_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qr_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qs_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qg_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_cf_it =>NULL() type(sub2grid_info) :: grd_t logical regional logical:: l_cld_derived,zflag,inithead type(gsi_bundle) :: atm_bundle type(gsi_grid) :: atm_grid integer(i_kind),parameter :: n2d=2 integer(i_kind),parameter :: n3d=14 character(len=4), parameter :: vars2d(n2d) = (/ 'z ', 'ps ' /) character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', & 'vor ', 'div ', & 'tv ', 'q ', & 'cw ', 'oz ', & 'ql ', 'qi ', & 'qr ', 'qs ', & 'qg ', 'cf ' /) real(r_kind),pointer,dimension(:,:):: ptr2d =>NULL() real(r_kind),pointer,dimension(:,:,:):: ptr3d =>NULL() regional=.false. inner_vars=1 ! num_fields=min(8*grd_a%nsig+2,npe) num_fields=min(n3d*grd_a%nsig+n2d,npe) if (mype==0) write(6,*)'npe num_fields = ', npe,n3d*grd_a%nsig+n2d,num_fields ! Create temporary communication information fore read routines call general_sub2grid_create_info(grd_t,inner_vars,grd_a%nlat,grd_a%nlon, & grd_a%nsig,num_fields,regional) ! Allocate bundle used for reading members call gsi_gridcreate(atm_grid,lat2,lon2,nsig) call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d) if(istatus/=0) then write(6,*) myname_,': trouble creating atm_bundle' call stop2(999) endif do it=1,nfldsig write(filename,'(''sigf'',i2.2)') ifilesig(it) ! Read background fields into bundle if (mype==0) write(6,*)'fv3_full_hydro = ', fv3_full_hydro if (fv3_full_hydro) then call general_read_fv3atm_nems(grd_t,sp_a,filename,.true.,.true.,.true.,& atm_bundle,.true.,istatus) else call general_read_gfsatm_nems(grd_t,sp_a,filename,.true.,.true.,.true.,& atm_bundle,.true.,istatus,it) endif inithead=.false. zflag=.false. ! Set values to actual MetGuess fields call set_guess_ if (it==ntguessig) then if (mype==0) write(6,*)'Print guess field ... after set_guess' call prt_guess('guess') endif if (fv3_full_hydro) then do k=1, nsig do j=1, lon2 do i=1, lat2 ! set lower bound to hydrometeors if (associated(ges_ql_it)) ges_ql_it(i,j,k) = max(qcmin,ges_ql_it(i,j,k)) if (associated(ges_qi_it)) ges_qi_it(i,j,k) = max(qcmin,ges_qi_it(i,j,k)) if (associated(ges_qr_it)) ges_qr_it(i,j,k) = max(qcmin,ges_qr_it(i,j,k)) if (associated(ges_qs_it)) ges_qs_it(i,j,k) = max(qcmin,ges_qs_it(i,j,k)) if (associated(ges_qg_it)) ges_qg_it(i,j,k) = max(qcmin,ges_qg_it(i,j,k)) if (associated(ges_cf_it)) ges_cf_it(i,j,k) = min(max(zero,ges_cf_it(i,j,k)),one) enddo enddo enddo else l_cld_derived = associated(ges_cwmr_it).and.& associated(ges_q_it) .and.& associated(ges_ql_it) .and.& associated(ges_qi_it) .and.& associated(ges_tv_it) ! call set_cloud_lower_bound(ges_cwmr_it) if (mype==0) write(6,*)'READ_GFS_NEMS: l_cld_derived = ', l_cld_derived if (l_cld_derived) then if (associated(ges_cf_it)) then call cloud_calc_gfs(ges_ql_it,ges_qi_it,ges_cwmr_it,ges_q_it,ges_tv_it,.true.,ges_cf_it) else call cloud_calc_gfs(ges_ql_it,ges_qi_it,ges_cwmr_it,ges_q_it,ges_tv_it,.true.) endif end if end if if (it==ntguessig) then if (mype==0) write(6,*)'Print guess field ... after reset cloud lower bound' call prt_guess('guess') endif end do call general_sub2grid_destroy_info(grd_t) call gsi_bundledestroy(atm_bundle,istatus) contains subroutine set_guess_ call gsi_bundlegetpointer (atm_bundle,'ps',ptr2d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ps',ges_ps_it ,istatus) if(istatus==0) ges_ps_it = ptr2d endif call gsi_bundlegetpointer (atm_bundle,'z',ptr2d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it ,istatus) if(istatus==0) ges_z_it = ptr2d endif call gsi_bundlegetpointer (atm_bundle,'u',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'u' ,ges_u_it ,istatus) if(istatus==0) ges_u_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'v',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'v' ,ges_v_it ,istatus) if(istatus==0) ges_v_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'vor',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'vor',ges_vor_it,istatus) if(istatus==0) ges_vor_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'div',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'div',ges_div_it,istatus) if(istatus==0) ges_div_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'tv',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'tv',ges_tv_it ,istatus) if(istatus==0) ges_tv_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'q',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'q' ,ges_q_it ,istatus) if(istatus==0) ges_q_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'oz',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'oz',ges_oz_it ,istatus) if(istatus==0) ges_oz_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'cw',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'cw',ges_cwmr_it,istatus) if(istatus==0) ges_cwmr_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'ql',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ql',ges_ql_it,istatus) if(istatus==0) ges_ql_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'qi',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qi',ges_qi_it,istatus) if(istatus==0) ges_qi_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'qr',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qr',ges_qr_it,istatus) if(istatus==0) ges_qr_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'qs',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qs',ges_qs_it,istatus) if(istatus==0) ges_qs_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'qg',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qg',ges_qg_it,istatus) if(istatus==0) ges_qg_it = ptr3d endif call gsi_bundlegetpointer (atm_bundle,'cf',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'cf',ges_cf_it,istatus) if(istatus==0) ges_cf_it = ptr3d endif end subroutine set_guess_ end subroutine read_ subroutine read_chem_ ( iyear, month,idd ) !$$$ subprogram documentation block ! . . . ! subprogram: read_nems_chem ! ! prgrmmr: todling ! ! abstract: fills chemguess_bundle with GFS chemistry. ! ! remarks: ! 1. Right now, only CO2 is done and even this is treated ! as constant througout the assimialation window. ! 2. iyear and month could come from obsmod, but logically ! this program should never depend on obsmod ! ! ! program history log: ! 2010-12-23 Huang - initial code, based on read_gfs_chem ! 2011-06-29 todling - no explict reference to internal bundle arrays ! 2019-04-19 Wei/Martin - modified to read NEMS aerosols from either ! NGAC or FV3-Chem ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: i_kind, r_kind use mpimod, only: mype use gridmod, only: lat2,lon2,nsig,nlat,rlats,istart use ncepgfs_ghg, only: read_gfsco2 use guess_grids, only: nfldsig use gsi_bundlemod, only: gsi_bundlegetpointer use gsi_chemguess_mod, only: gsi_chemguess_bundle use gsi_chemguess_mod, only: gsi_chemguess_get use gsi_bundlemod, only: gsi_bundle,gsi_bundlecreate,gsi_bundledestroy use gsi_bundlemod, only: gsi_grid,gsi_gridcreate use radiance_mod, only: n_aerosols_fwd,aerosol_names_fwd use gridmod, only: grd_a,sp_a,regional use guess_grids, only: ifilesig,ifileaer,nfldaer use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info use mpimod, only: npe use chemmod, only: lread_ext_aerosol implicit none ! Declared argument list integer(i_kind), intent(in):: iyear integer(i_kind), intent(in):: month integer(i_kind), intent(in):: idd ! Declare local variables integer(i_kind) :: igfsco2, i, j, n, iret real(r_kind),dimension(lat2):: xlats real(r_kind),pointer,dimension(:,:,:)::p_co2=>NULL() real(r_kind),pointer,dimension(:,:,:)::ptr3d=>NULL() integer(i_kind) :: i4crtm real(r_kind),pointer,dimension(:,:,:)::ae_du001_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_du002_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_du003_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_du004_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_du005_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_ss001_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_ss002_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_ss003_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_ss004_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_so4_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_ocpho_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_ocphi_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_bcpho_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ae_bcphi_it=>NULL() type(sub2grid_info) :: grd_ae character(24) filename,str_crtmuse integer(i_kind) :: it,ia,ier,iera,inner_vars,num_fields,istatus integer(i_kind) :: anfld type(gsi_bundle) :: chem_bundle type(gsi_grid) :: chem_grid if(.not.associated(gsi_chemguess_bundle)) return call gsi_bundlegetpointer(gsi_chemguess_bundle(1),'co2',p_co2,iret) if(iret /= 0) return ! Get subdomain latitude array j = mype + 1 do i = 1, lat2 n = min(max(1, istart(j)+i-2), nlat) xlats(i) = rlats(n) enddo ! Read in CO2 call gsi_chemguess_get ( 'i4crtm::co2', igfsco2, iret ) call read_gfsco2 ( iyear,month,idd,igfsco2,xlats,& lat2,lon2,nsig,mype, p_co2 ) ! Approximation: setting all times co2 values equal to the daily co2 values do n = 2, nfldsig call gsi_bundlegetpointer(gsi_chemguess_bundle(n),'co2',ptr3d,iret) ptr3d = p_co2 enddo ! Read in Aerosol field via nemsio if ( n_aerosols_fwd > 0 ) then if ( mype == 0 ) write(6,*) 'n_aerosols_fwd and aerosol_names_fwd',n_aerosols_fwd,aerosol_names_fwd call gsi_gridcreate(chem_grid,lat2,lon2,nsig) call gsi_bundlecreate(chem_bundle,chem_grid,'aux-chem-read',istatus,names3d=aerosol_names_fwd) inner_vars=1 num_fields=min(n_aerosols_fwd*grd_a%nsig,npe) ! Create temporary communication information fore read routines call general_sub2grid_create_info(grd_ae,inner_vars,grd_a%nlat,grd_a%nlon,grd_a%nsig,num_fields,regional) if (lread_ext_aerosol) then anfld=nfldaer else anfld=nfldsig end if do it=1,anfld ! Get pointer to aerosol field if (lread_ext_aerosol) then write(filename,'(''aerf'',i2.2)') ifileaer(it) else write(filename,'(''sigf'',i2.2)') ifilesig(it) end if if (mype==0) write(6,*) "aerosol field come from ",filename ier=0 call general_read_nemsaero(grd_ae,sp_a,filename,mype,chem_bundle,& n_aerosols_fwd,aerosol_names_fwd,.true.,ier) do ia=1,n_aerosols_fwd write(str_crtmuse,'(''i4crtm::'',a)') trim(aerosol_names_fwd(ia)) call gsi_chemguess_get ( str_crtmuse, i4crtm, iera ) if(mype==0) write(6,*) trim(aerosol_names_fwd(ia))," for crtm is ",i4crtm call gsi_bundlegetpointer (chem_bundle,trim(aerosol_names_fwd(ia)),ptr3d,istatus) if (istatus==0) then select case ( trim(aerosol_names_fwd(ia)) ) case ('sulf') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'sulf',ae_so4_it,iret) if (iret==0) ae_so4_it=ptr3d case ('bc1') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'bc1',ae_bcpho_it,iret) if (iret==0) ae_bcpho_it=ptr3d case ('bc2') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'bc2',ae_bcphi_it,iret) if (iret==0) ae_bcphi_it=ptr3d case ('oc1') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'oc1',ae_ocpho_it,iret) if (iret==0) ae_ocpho_it=ptr3d case ('oc2') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'oc2',ae_ocphi_it,iret) if (iret==0) ae_ocphi_it=ptr3d case ('dust1') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'dust1',ae_du001_it,iret) if (iret==0) ae_du001_it=ptr3d case ('dust2') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'dust2',ae_du002_it,iret) if (iret==0) ae_du002_it=ptr3d case ('dust3') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'dust3',ae_du003_it,iret) if (iret==0) ae_du003_it=ptr3d case ('dust4') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'dust4',ae_du004_it,iret) if (iret==0) ae_du004_it=ptr3d case ('dust5') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'dust5',ae_du005_it,iret) if (iret==0) ae_du005_it=ptr3d case ('seas1') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'seas1',ae_ss001_it,iret) if (iret==0) ae_ss001_it=ptr3d case ('seas2') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'seas2',ae_ss002_it,iret) if (iret==0) ae_ss002_it=ptr3d case ('seas3') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'seas3',ae_ss003_it,iret) if (iret==0) ae_ss003_it=ptr3d case ('seas4') call gsi_bundlegetpointer(gsi_chemguess_bundle(it),'seas4',ae_ss004_it,iret) if (iret==0) ae_ss004_it=ptr3d end select ! different aerosol tracers if(iret/=0 .and. mype==0 ) write(6,*) trim(aerosol_names_fwd(ia))," getpointer fail" ier=ier+iret endif ! end if successfully able to get pointer end do ! n_aerosols_fwd if (ier/=0) then write(6,*) "before call read_ngac_aerosol ier=",ier cycle ! this allows code to be free from met-fields end if end do ! nfldaer call general_sub2grid_destroy_info(grd_ae) call gsi_bundledestroy(chem_bundle,istatus) end if ! end if n_aerosols_fwd > 0 end subroutine read_chem_ subroutine read_atm_ (grd,filename,sp_a,uvflag,vordivflag,zflag, & g_z,g_ps,g_vor,g_div,g_u,g_v,& g_tv,g_q,g_cwmr,g_oz) !$$$ subprogram documentation block ! . . . . ! subprogram: read_nemsatm read nems atm and send to all mpi tasks ! prgmmr: Huang org: np23 date: 2010-02-22 ! ! abstract: read ncep nems/gfs atmospheric guess field and ! scatter to subdomains ! ! program history log: ! 2010-02-22 Huang Initial version. Based on sub read_gfsatm ! 2011-02-28 Huang Re-arrange the read sequence to be same as model ! write sequence. Alsom allocate and deallocate ! temporary working array immediatelt before and after ! the processing and scattering first guess field to reduce ! maximum resident memory size. Page fault can happen ! when running at high resolution GSI, e.g., T574. ! 2011-09-23 Huang Add NEMS parallel IO capability ! 2013-10-25 todling reposition fill_ns,filluv_ns to commvars ! ! input argument list: ! grd - structure variable containing information about grid ! (initialized by general_sub2grid_create_info, located in ! general_sub2grid_mod.f90) ! sp_a - structure variable containing spectral information for analysis ! (initialized by general_init_spec_vars, located in ! general_specmod.f90) ! uvflag - logical to use u,v (.true.) or st,vp (.false.) perturbations ! vordivflag - logical to determine if routine should output vorticity and ! divergence ! zflag - logical to determine if surface height field should be output ! ! output argument list: ! g_* - guess fields ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind use gridmod, only: ntracer,ncloud,reload,itotsub use general_commvars_mod, only: fill_ns,filluv_ns,fill2_ns,filluv2_ns,ltosj_s,ltosi_s use general_specmod, only: spec_vars use general_sub2grid_mod, only: sub2grid_info use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype,mype use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close use nemsio_module, only: nemsio_gfile,nemsio_getfilehead,nemsio_readrecv use egrid2agrid_mod,only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid use constants, only: two,pi,half,deg2rad implicit none ! Declare local parameters real(r_kind),parameter:: r0_001 = 0.001_r_kind ! Declare passed variables type(sub2grid_info) ,intent(in ) :: grd character(len=24) ,intent(in ) :: filename logical ,intent(in ) :: uvflag,vordivflag,zflag real(r_kind),dimension(grd%lat2,grd%lon2) ,intent( out) :: g_z,g_ps real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( out) :: g_u,g_v,& g_vor,g_div,g_cwmr,g_q,g_oz,g_tv type(spec_vars) ,intent(in ) :: sp_a ! Declare local variables character(len=120) :: my_name = 'READ_NEMSATM' character(len=1) :: null = ' ' integer(i_kind),dimension(7):: idate integer(i_kind),dimension(4):: odate integer(i_kind) :: iret,nlatm2,nflds integer(i_kind) :: k,icount,icount_prev,mm1,i,j,kk integer(i_kind) :: mype_hs, mype_ps,nord_int integer(i_kind) :: latb, lonb, levs, nframe integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd integer(i_kind) :: istop = 101 real(r_kind),allocatable,dimension(:,:) :: grid, grid_v, & grid_vor, grid_div, grid_b, grid_b2 real(r_kind),allocatable,dimension(:,:,:) :: grid_c, grid2, grid_c2 real(r_kind),allocatable,dimension(:) :: work, work_vor, work_div, & work_v real(r_kind),allocatable,dimension(:,:) :: sub, sub_vor, sub_div, & sub_v real(r_kind),dimension(sp_a%nc):: spec_vor,spec_div real(r_kind),allocatable,dimension(:) :: rwork1d0, rwork1d1, rwork1d2 real(r_kind),allocatable,dimension(:) :: rlats,rlons,clons,slons real(4),allocatable,dimension(:) :: r4lats,r4lons real(r_kind) :: fhour type(nemsio_gfile) :: gfile logical diff_res,eqspace logical,dimension(1) :: vector type(egrid2agrid_parm) :: p_high !****************************************************************************** ! Initialize variables used below mm1=mype+1 mype_hs=min(1,npe-1) mype_ps=0 nlatm2=grd%nlat-2 nflds=5*grd%nsig+1 if(zflag) nflds=nflds+1 if(vordivflag .or. .not. uvflag)nflds=nflds+2*grd%nsig ! nflds=npe nflds=grd%nsig levs=grd%nsig allocate( work(grd%itotsub),work_v(grd%itotsub) ) work=zero work_v=zero allocate( sub(grd%lat2*grd%lon2,max(grd%nsig,npe)),sub_v(grd%lat2*grd%lon2,max(grd%nsig,npe)) ) allocate( sub_div(grd%lat2*grd%lon2,max(grd%nsig,npe)),sub_vor(grd%lat2*grd%lon2,max(grd%nsig,npe)) ) if(mype < nflds)then call nemsio_init(iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'init',istop,iret) call nemsio_open(gfile,filename,'READ',iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop+1,iret) call nemsio_getfilehead(gfile,iret=iret, nframe=nframe, & nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, & idate=idate, dimx=lonb, dimy=latb,dimz=levs) if( nframe /= 0 ) then if ( mype == 0 ) & write(6,*)trim(my_name),': ***ERROR*** nframe /= 0 for global model read, nframe = ', nframe call stop2(101) end if fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600 odate(1) = idate(4) !hour odate(2) = idate(2) !month odate(3) = idate(3) !day odate(4) = idate(1) !year ! ! g_* array already pre-allocate as (lat2,lon2,) => 2D and <3D> array ! diff_res=.false. if(latb /= nlatm2) then diff_res=.true. if ( mype == 0 ) write(6, & '(a,'': different spatial dimension nlatm2 = '',i4,tr1,''latb = '',i4)') & trim(my_name),nlatm2,latb ! call stop2(101) end if if(lonb /= grd%nlon) then diff_res=.true. if ( mype == 0 ) write(6, & '(a,'': different spatial dimension nlon = '',i4,tr1,''lonb = '',i4)') & trim(my_name),grd%nlon,lonb ! call stop2(101) end if if(levs /= grd%nsig)then if ( mype == 0 ) write(6, & '(a,'': inconsistent spatial dimension nsig = '',i4,tr1,''levs = '',i4)') & trim(my_name),grd%nsig,levs call stop2(101) end if ! allocate( grid(grd%nlon,nlatm2), grid_v(grd%nlon,nlatm2) ) if(diff_res)then allocate( grid_b(lonb,latb),grid_c(latb+2,lonb,1),grid2(grd%nlat,grd%nlon,1)) allocate( grid_b2(lonb,latb),grid_c2(latb+2,lonb,1)) end if allocate( rwork1d0(latb*lonb) ) allocate( rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb),r4lats(lonb*latb),r4lons(lonb*latb)) allocate(rwork1d1(latb*lonb)) call nemsio_getfilehead(gfile,lat=r4lats,iret=iret) call nemsio_getfilehead(gfile,lon=r4lons,iret=iret) do j=1,latb rlats(latb+2-j)=deg2rad*r4lats(lonb/2+(j-1)*lonb) end do do j=1,lonb rlons(j)=deg2rad*r4lons(j) end do deallocate(r4lats,r4lons) rlats(1)=-half*pi rlats(latb+2)=half*pi do j=1,lonb clons(j)=cos(rlons(j)) slons(j)=sin(rlons(j)) end do nord_int=4 eqspace=.false. call g_create_egrid2agrid(grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons, & latb+2,rlats,lonb,rlons,& nord_int,p_high,.true.,eqspace) deallocate(rlats,rlons) end if ! ! Load values into rows for south and north pole before scattering ! ! Terrain: scatter to all mpi tasks ! if(zflag)then if (mype==mype_hs) then call nemsio_readrecv(gfile,'hgt', 'sfc',1,rwork1d0,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'hgt','read',istop+2,iret) if(diff_res)then grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) do kk=1,itotsub i=ltosi_s(kk) j=ltosj_s(kk) work(kk)=grid2(i,j,1) end do else grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/)) call fill_ns(grid,work) end if endif call mpi_scatterv(work,grd%ijn_s,grd%displs_s,mpi_rtype,& g_z,grd%ijn_s(mm1),mpi_rtype,mype_hs,mpi_comm_world,ierror) end if ! Surface pressure: same procedure as terrain, but handled by task mype_ps ! if (mype==mype_ps) then call nemsio_readrecv(gfile,'pres','sfc',1,rwork1d0,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pres','read',istop+3,iret) rwork1d1 = r0_001*rwork1d0 if(diff_res)then vector(1)=.false. grid_b=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/)) call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) do kk=1,itotsub i=ltosi_s(kk) j=ltosj_s(kk) work(kk)=grid2(i,j,1) end do else grid=reshape(rwork1d1,(/size(grid,1),size(grid,2)/)) ! convert Pa to cb call fill_ns(grid,work) endif endif call mpi_scatterv(work,grd%ijn_s,grd%displs_s,mpi_rtype,& g_ps,grd%ijn_s(mm1),mpi_rtype,mype_ps,mpi_comm_world,ierror) ! Divergence and voriticity. Compute u and v from div and vor sub_vor=zero sub_div=zero sub =zero sub_v =zero icount =0 icount_prev=1 allocate( work_vor(grd%itotsub),work_div(grd%itotsub) ) do k=1,levs icount=icount+1 if (mype==mod(icount-1,npe)) then ! Convert grid u,v to div and vor call nemsio_readrecv(gfile,'ugrd','mid layer',k,rwork1d0,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ugrd','read',istop+4,iret) call nemsio_readrecv(gfile,'vgrd','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vgrd','read',istop+5,iret) if(diff_res)then grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/)) grid_b2=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.true. call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) do kk=1,itotsub i=ltosi_s(kk) j=ltosj_s(kk) work(kk)=grid2(i,j,1) end do do j=1,grd%nlon do i=2,grd%nlat-1 grid(j,grd%nlat-i)=grid2(i,j,1) end do end do call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) do kk=1,itotsub i=ltosi_s(kk) j=ltosj_s(kk) work_v(kk)=grid2(i,j,1) end do do j=1,grd%nlon do i=2,grd%nlat-1 grid_v(j,grd%nlat-i)=grid2(i,j,1) end do end do else grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/)) grid_v=reshape(rwork1d1,(/size(grid_v,1),size(grid_v,2)/)) call filluv_ns(grid,grid_v,work,work_v) end if if(vordivflag .or. .not. uvflag)then allocate( grid_vor(grd%nlon,nlatm2), grid_div(grd%nlon,nlatm2) ) call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) call general_sptez_s_b(sp_a,sp_a,spec_div,grid_div,1) call general_sptez_s_b(sp_a,sp_a,spec_vor,grid_vor,1) ! Load values into rows for south and north pole call fill_ns(grid_div,work_div) call fill_ns(grid_vor,work_vor) deallocate(grid_vor,grid_div) end if endif ! Scatter to sub if (mod(icount,npe)==0 .or. icount==levs) then if(vordivflag .or. .not. uvflag)then call mpi_alltoallv(work_vor,grd%ijn_s,grd%displs_s,mpi_rtype,& sub_vor(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,& mpi_comm_world,ierror) call mpi_alltoallv(work_div,grd%ijn_s,grd%displs_s,mpi_rtype,& sub_div(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,& mpi_comm_world,ierror) end if if(uvflag)then call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,& sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,& mpi_comm_world,ierror) call mpi_alltoallv(work_v,grd%ijn_s,grd%displs_s,mpi_rtype,& sub_v(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,& mpi_comm_world,ierror) end if icount_prev=icount+1 endif end do deallocate(work_vor,work_div) ! Transfer vor,div,u,v into real(r_kind) guess arrays call reload(sub_vor,g_vor) call reload(sub_div,g_div) call reload(sub,g_u) call reload(sub_v,g_v) deallocate(sub_vor,sub_div) ! Thermodynamic variable and Specific humidity: communicate to all tasks ! sub=zero icount=0 icount_prev=1 do k=1,levs icount=icount+1 if (mype==mod(icount-1,npe)) then call nemsio_readrecv(gfile,'spfh','mid layer',k,rwork1d0,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','read',istop+6,iret) if(diff_res)then grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) do kk=1,itotsub i=ltosi_s(kk) j=ltosj_s(kk) work(kk)=grid2(i,j,1) end do else grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/)) call fill_ns(grid,work) end if call nemsio_readrecv(gfile,'tmp','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','read',istop+7,iret) allocate(rwork1d2(latb*lonb)) rwork1d2 = rwork1d1*(one+fv*rwork1d0) if(diff_res)then grid_b=reshape(rwork1d2,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) do kk=1,itotsub i=ltosi_s(kk) j=ltosj_s(kk) work_v(kk)=grid2(i,j,1) end do else grid_v=reshape(rwork1d2,(/size(grid_v,1),size(grid_v,2)/)) call fill_ns(grid_v,work_v) end if deallocate(rwork1d2) endif if (mod(icount,npe)==0 .or. icount==levs) then call mpi_alltoallv(work_v,grd%ijn_s,grd%displs_s,mpi_rtype,& sub_v(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,& mpi_comm_world,ierror) call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,& sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,& mpi_comm_world,ierror) icount_prev=icount+1 endif end do call reload(sub_v,g_tv) call reload(sub,g_q) deallocate(sub_v,work_v) sub=zero icount=0 icount_prev=1 do k=1,levs icount=icount+1 if (mype==mod(icount-1,npe)) then call nemsio_readrecv(gfile,'o3mr','mid layer',k,rwork1d0,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','read',istop+8,iret) if(diff_res)then grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) do kk=1,itotsub i=ltosi_s(kk) j=ltosj_s(kk) work(kk)=grid2(i,j,1) end do else grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/)) call fill_ns(grid,work) end if endif if (mod(icount,npe)==0 .or. icount==levs) then call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,& sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,& mpi_comm_world,ierror) icount_prev=icount+1 endif end do call reload(sub,g_oz) ! Cloud condensate mixing ratio. if (ntracer>2 .or. ncloud>=1) then sub=zero icount=0 icount_prev=1 do k=1,levs icount=icount+1 if (mype==mod(icount-1,npe)) then call nemsio_readrecv(gfile,'clwmr','mid layer',k,rwork1d0,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','read',istop+9,iret) if (imp_physics == 11) then call nemsio_readrecv(gfile,'icmr','mid layer',k,rwork1d1,iret=iret) if (iret == 0) then rwork1d0 = rwork1d0 + rwork1d1 else call error_msg(trim(my_name),trim(filename),'icmr','read',istop+10,iret) endif endif if(diff_res)then grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) do kk=1,itotsub i=ltosi_s(kk) j=ltosj_s(kk) work(kk)=grid2(i,j,1) end do else grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/)) call fill_ns(grid,work) end if endif if (mod(icount,npe)==0 .or. icount==levs) then call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,& sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,& mpi_comm_world,ierror) icount_prev=icount+1 endif end do call reload(sub,g_cwmr) else g_cwmr = zero endif if(mype < nflds)then if(diff_res) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid2) call destroy_egrid2agrid(p_high) deallocate(rwork1d1,clons,slons) deallocate(rwork1d0) deallocate(grid,grid_v) call nemsio_close(gfile,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop+9,iret) end if deallocate(work,sub) ! Print date/time stamp if ( mype == 0 ) write(6, & '(a,'': ges read/scatter,lonb,latb,levs= '',3i6,'',hour= '',f4.1,'',idate= '',4i5)') & trim(my_name),lonb,latb,levs,fhour,odate end subroutine read_atm_ subroutine read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, & veg_type,soil_type,terrain,isli,use_sfc_any, & tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d) !$$$ subprogram documentation block ! . . . . ! subprogram: read_sfc_ read nems hist file ! prgmmr: Li org: np23 date: 2017-08-30 ! ! abstract: read nems sfc & nst combined file ! ! program history log: ! ! input argument list: ! use_sfc_any - true if any processor uses extra surface fields ! ! output argument list: ! sfct - surface temperature (skin temp) ! soil_moi - soil moisture of first layer ! sno - snow depth ! soil_temp - soil temperature of first layer ! veg_frac - vegetation fraction ! fact10 - 10 meter wind factor ! sfc_rough - surface roughness ! veg_type - vegetation type ! soil_type - soil type ! terrain - terrain height ! isli - sea/land/ice mask ! tref - optional, oceanic foundation temperature ! dt_cool - optional, sub-layer cooling amount at sub-skin layer ! z_c - optional, depth of sub-layer cooling layer ! dt_warm - optional, diurnal warming amount at sea surface ! z_w - optional, depth of diurnal warming layer ! c_0 - optional, coefficient to calculate d(Tz)/d(tr) in dimensionless ! c_d - optional, coefficient to calculate d(Tz)/d(tr) in m^-1 ! w_0 - optional, coefficient to calculate d(Tz)/d(tr) in dimensionless ! w_d - optional, coefficient to calculate d(Tz)/d(tr) in m^-1 ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use mpimod, only: mype use kinds, only: r_kind,i_kind,r_single use gridmod, only: nlat_sfc,nlon_sfc use guess_grids, only: nfldsfc,ifilesfc use constants, only: zero,two use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close use nemsio_module, only: nemsio_gfile,nemsio_getfilehead,nemsio_readrecv implicit none ! Declare passed variables logical, intent(in) :: use_sfc_any real(r_single), dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough real(r_single), dimension(nlat_sfc,nlon_sfc), intent(out) :: veg_type,soil_type,terrain integer(i_kind), dimension(nlat_sfc,nlon_sfc), intent(out) :: isli real(r_single), optional, dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d ! Declare local parameters integer(i_kind), parameter :: nsfc_all=11 integer(i_kind),dimension(7):: idate integer(i_kind),dimension(4):: odate ! Declare local variables real(r_single), dimension(nlat_sfc,nlon_sfc,nfldsfc) :: xt character(len=24) :: filename character(len=120) :: my_name = 'READ_SFCNST' character(len=1) :: null = ' ' integer(i_kind) :: i,j,it,n,nsfc integer(i_kind) :: iret, nframe, lonb, latb integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd real(r_single) :: fhour integer(i_kind) :: istop = 102 real(r_single), allocatable, dimension(:) :: rwork2d real(r_single), allocatable, dimension(:,:) :: work,outtmp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Define read variable property !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! type(nemsio_gfile) :: gfile !----------------------------------------------------------------------------- call nemsio_init(iret=iret) if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret) do it = 1, nfldsfc ! read a surface history file on the task write(filename,200)ifilesfc(it) 200 format('sfcf',i2.2) call nemsio_open(gfile,filename,'READ',iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret) call nemsio_getfilehead(gfile, idate=idate, iret=iret, nframe=nframe, & nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, & dimx=lonb, dimy=latb ) if( nframe /= 0 ) then if ( mype == 0 ) & write(6,*)trim(my_name),': ***ERROR*** nframe /= 0 for global sfc hist read, nframe = ', nframe call stop2(102) end if fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600 odate(1) = idate(4) !hour odate(2) = idate(2) !month odate(3) = idate(3) !day odate(4) = idate(1) !year if ( (latb /= nlat_sfc-2) .or. (lonb /= nlon_sfc) ) then if ( mype == 0 ) write(6, & '(a,'': inconsistent spatial dimension '',''nlon,nlatm2 = '',2(i4,tr1),''-vs- sfc file lonb,latb = '',i4)') & trim(my_name),nlon_sfc,nlat_sfc-2,lonb,latb call stop2(102) endif ! ! Read the surface records (lonb, latb) and convert to GSI array pattern (nlat_sfc,nlon_sfc) ! Follow the read order sfcio in ncepgfs_io ! allocate(work(lonb,latb)) allocate(rwork2d(size(work,1)*size(work,2))) work = zero rwork2d = zero if(it == 1)then nsfc=nsfc_all else nsfc=nsfc_all-4 end if do n = 1, nsfc if (n == 1) then ! skin temperature ! Tsea call nemsio_readrecv(gfile, 'tmp', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,sfct(1,1,it),lonb,latb) elseif(n == 2 .and. use_sfc_any) then ! soil moisture ! smc/soilw call nemsio_readrecv(gfile, 'smc', 'soil layer', 1, rwork2d, iret=iret) ! FV3 nemsio files use 'soilw 0-10cm down' insted of 'smc soil layer 1' if (iret /= 0) then if ( mype == 0 ) print *,'could not read smc, try to read soilw 0-10 cm down instead...' call nemsio_readrecv(gfile,'soilw','0-10 cm down',1,rwork2d,iret=iret) if (iret /= 0) & call error_msg(trim(my_name),trim(filename),'smc/soilw','read',istop,iret) endif work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,soil_moi(1,1,it),lonb,latb) elseif(n == 3) then ! snow depth ! sheleg call nemsio_readrecv(gfile, 'weasd','sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'weasd','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,sno(1,1,it),lonb,latb) elseif(n == 4 .and. use_sfc_any) then ! soil temperature ! stc/tmp call nemsio_readrecv(gfile, 'stc', 'soil layer', 1, rwork2d, iret=iret) if (iret /= 0) then ! FV3 nemsio files use 'tmp 0-10cm down' insted of 'stc soil layer 1' if ( mype == 0 ) print *,'could not read stc, try to read tmp 0-10 cm down instead...' call nemsio_readrecv(gfile,'tmp','0-10 cm down',1,rwork2d,iret=iret) if (iret /= 0) & call error_msg(trim(my_name),trim(filename),'stc/tmp','read',istop,iret) endif work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,soil_temp(1,1,it),lonb,latb) elseif(n == 5 .and. use_sfc_any) then ! vegetation cover ! vfrac call nemsio_readrecv(gfile, 'veg', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'veg','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,veg_frac(1,1,it),lonb,latb) elseif(n == 6) then ! 10m wind factor ! f10m call nemsio_readrecv(gfile, 'f10m', '10 m above gnd', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'f10m','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,fact10(1,1,it),lonb,latb) elseif(n == 7) then ! suface roughness ! zorl call nemsio_readrecv(gfile, 'sfcr', 'sfc', 1, rwork2d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'sfcr','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,sfc_rough(1,1,it),lonb,latb) elseif(n == 8 .and. use_sfc_any) then ! vegetation type ! vtype call nemsio_readrecv(gfile, 'vtype','sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vtype','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,veg_type,lonb,latb) elseif(n == 9 .and. use_sfc_any) then ! soil type ! stype call nemsio_readrecv(gfile, 'sotyp','sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'sotyp','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,soil_type,lonb,latb) elseif(n == 10) then ! terrain ! orog call nemsio_readrecv(gfile, 'orog', 'sfc', 1, rwork2d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'orog','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,terrain,lonb,latb) elseif(n == 11) then ! sea/land/ice flag ! slmsk call nemsio_readrecv(gfile, 'land', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'land','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) allocate(outtmp(latb+2,lonb)) call tran_gfssfc(work,outtmp,lonb,latb) do j=1,lonb do i=1,latb+2 isli(i,j) = nint(outtmp(i,j)) end do end do deallocate(outtmp) endif ! End of loop over data records enddo if( present(tref) ) then if ( mype == 0 ) write(6,*) ' read 9 optional NSST variables ' call nemsio_readrecv(gfile, 'tref', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tref','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,tref(1,1,it),lonb,latb) call nemsio_readrecv(gfile, 'dtcool','sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dtcool','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,dt_cool(1,1,it),lonb,latb) call nemsio_readrecv(gfile, 'zc', 'sfc', 1, rwork2d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'z_c','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,z_c(1,1,it),lonb,latb) call nemsio_readrecv(gfile, 'xt', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'xt','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,xt(1,1,it),lonb,latb) call nemsio_readrecv(gfile, 'xz', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'xz','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,z_w(1,1,it),lonb,latb) call nemsio_readrecv(gfile, 'c0', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'c0','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,c_0(1,1,it),lonb,latb) call nemsio_readrecv(gfile, 'cd', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'cd','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,c_d(1,1,it),lonb,latb) call nemsio_readrecv(gfile, 'w0', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'w0','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,w_0(1,1,it),lonb,latb) call nemsio_readrecv(gfile, 'wd', 'sfc', 1, rwork2d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'wd','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,w_d(1,1,it),lonb,latb) ! ! Get diurnal warming amout at z=0 ! do j = 1,nlon_sfc do i = 1,nlat_sfc if (z_w(i,j,it) > zero) then dt_warm(i,j,it) = two*xt(i,j,it)/z_w(i,j,it) end if end do end do endif ! Deallocate local work arrays deallocate(work,rwork2d) call nemsio_close(gfile,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret) ! ! Print date/time stamp if ( mype == 0 ) write(6, & '(a,'': sfc read,nlon,nlat= '',2i6,'',hour= '',f4.1,'',idate= '',4i5)') & trim(my_name),lonb,latb,fhour,odate ! End of loop over time levels end do end subroutine read_sfc_ subroutine read_nemssfc_(iope,sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, & veg_type,soil_type,terrain,isli,use_sfc_any, & tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d) !$$$ subprogram documentation block ! . . . . ! subprogram: read_nemssfc_ read nems hist file ! prgmmr: xuli org: np23 date: 2017-08-30 ! ! abstract: read nems surface file ! ! program history log: ! 2017-08-30 li ! ! input argument list: ! iope - mpi task handling i/o ! use_sfc_any - true if any processor uses extra surface fields ! ! output argument list: ! sfct - surface temperature (skin temp) ! soil_moi - soil moisture of first layer ! sno - snow depth ! soil_temp - soil temperature of first layer ! veg_frac - vegetation fraction ! fact10 - 10 meter wind factor ! sfc_rough - surface roughness ! veg_type - vegetation type ! soil_type - soil type ! terrain - terrain height ! isli - sea/land/ice mask ! tref - oceanic foundation temperature ! dt_cool - optional, sub-layer cooling amount at sub-skin layer ! z_c - optional, depth of sub-layer cooling layer ! dt_warm - optional, diurnal warming amount at sea surface ! z_w - optional, depth of diurnal warming layer ! c_0 - optional, coefficient to calculate d(Tz)/d(tf) ! c_d - optional, coefficient to calculate d(Tz)/d(tf) ! w_0 - optional, coefficient to calculate d(Tz)/d(tf) ! w_d - optional, coefficient to calculate d(Tz)/d(tf) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind,r_single use gridmod, only: nlat_sfc,nlon_sfc use guess_grids, only: nfldsfc,sfcmod_mm5,sfcmod_gfs use mpimod, only: mpi_itype,mpi_rtype4,mpi_comm_world,mype use constants, only: zero implicit none ! Declare passed variables integer(i_kind), intent(in) :: iope logical, intent(in) :: use_sfc_any real(r_single), dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough real(r_single), dimension(nlat_sfc,nlon_sfc), intent(out) :: veg_type,soil_type,terrain integer(i_kind), dimension(nlat_sfc,nlon_sfc), intent(out) :: isli real(r_single), optional, dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d ! Declare local variables integer(i_kind):: iret,npts,nptsall !----------------------------------------------------------------------------- ! Read surface history file on processor iope if(mype == iope)then if ( present(tref) ) then call read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, & veg_type,soil_type,terrain,isli,use_sfc_any, & tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d) write(*,*) 'read_sfc nemsio, with NSST variables' else call read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, & veg_type,soil_type,terrain,isli,use_sfc_any) write(*,*) 'read_sfc nemsio, without NSST variables' endif endif ! Load onto all processors npts=nlat_sfc*nlon_sfc nptsall=npts*nfldsfc call mpi_bcast(sfct, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(fact10, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(sno, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) if(sfcmod_mm5 .or. sfcmod_gfs)then call mpi_bcast(sfc_rough, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) else sfc_rough = zero endif call mpi_bcast(terrain, npts, mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(isli, npts, mpi_itype, iope,mpi_comm_world,iret) if(use_sfc_any)then call mpi_bcast(veg_frac, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(soil_temp,nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(soil_moi, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(veg_type, npts, mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(soil_type,npts, mpi_rtype4,iope,mpi_comm_world,iret) endif if ( present(tref) ) then call mpi_bcast(tref, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(dt_cool, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(z_c, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(dt_warm, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(z_w, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(c_0, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(c_d, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(w_0, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(w_d, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) endif end subroutine read_nemssfc_ subroutine read_sfc_anl_(isli_anl) !$$$ subprogram documentation block ! . . . . ! subprogram: read_sfc_anl_ read nems surface file with analysis resolution ! ! prgmmr: li org: np23 date: 2016-08-18 ! ! abstract: read nems surface file at analysis grids when nlon /= nlon_sfc or nlat /= nlat_sfc ! ! program history log: ! ! input argument list: ! ! output argument list: ! isli - sea/land/ice mask ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use mpimod, only: mype use kinds, only: r_kind,i_kind,r_single use gridmod, only: nlat,nlon use constants, only: zero use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close use nemsio_module, only: nemsio_gfile,nemsio_getfilehead,nemsio_readrecv implicit none ! Declare passed variables integer(i_kind), dimension(nlat,nlon), intent( out) :: isli_anl ! Declare local parameters integer(i_kind),dimension(7):: idate integer(i_kind),dimension(4):: odate ! Declare local variables character(len=24) :: filename character(len=120) :: my_name = 'READ_NEMSSFC_ANL' character(len=1) :: null = ' ' integer(i_kind) :: i,j integer(i_kind) :: iret, nframe, lonb, latb integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd real(r_single) :: fhour integer(i_kind) :: istop = 102 real(r_single), allocatable, dimension(:) :: rwork2d real(r_single), allocatable, dimension(:,:) :: work,outtmp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Define read variable property !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! type(nemsio_gfile) :: gfile !----------------------------------------------------------------------------- call nemsio_init(iret=iret) if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret) filename='sfcf06_anlgrid' call nemsio_open(gfile,trim(filename),'READ',iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret) call nemsio_getfilehead(gfile, idate=idate, iret=iret, nframe=nframe, & nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, & dimx=lonb, dimy=latb ) if( nframe /= 0 ) then if ( mype == 0 ) & write(6,*)trim(my_name),': ***ERROR*** nframe /= 0 for global model read, nframe = ', nframe call stop2(102) end if fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600 odate(1) = idate(4) !hour odate(2) = idate(2) !month odate(3) = idate(3) !day odate(4) = idate(1) !year if ( (latb /= nlat-2) .or. (lonb /= nlon) ) then if ( mype == 0 ) write(6, & '(a,'': inconsistent spatial dimension '',''nlon,nlatm2 = '',2(i4,tr1),''-vs- sfc file lonb,latb = '',i4)') & trim(my_name),nlon,nlat-2,lonb,latb call stop2(102) endif ! ! Read the surface records (lonb, latb) and convert to GSI array pattern (nlat,nlon) ! Follow the read order sfcio in ncepgfs_io ! allocate(work(lonb,latb)) allocate(rwork2d(size(work,1)*size(work,2))) work = zero rwork2d = zero ! slmsk call nemsio_readrecv(gfile, 'land', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'land','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) allocate(outtmp(latb+2,lonb)) call tran_gfssfc(work,outtmp,lonb,latb) do j=1,lonb do i=1,latb+2 isli_anl(i,j) = nint(outtmp(i,j)) end do end do deallocate(outtmp) ! Deallocate local work arrays deallocate(work,rwork2d) call nemsio_close(gfile,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret) ! ! Print date/time stamp if ( mype == 0 ) write(6, & '(a,'': read_sfc_anl_ ,nlon,nlat= '',2i6,'',hour= '',f4.1,'',idate= '',4i5)') & trim(my_name),lonb,latb,fhour,odate end subroutine read_sfc_anl_ subroutine read_nemssfc_anl_(iope,isli_anl) !$$$ subprogram documentation block ! . . . . ! subprogram: read_nemssfc_anl read nems surface guess file with analysis resolution ! ! prgmmr: xuli org: np23 date: 2016-08-18 ! ! abstract: read nems surface file at analysis grids ! ! program history log: ! ! input argument list: ! iope - mpi task handling i/o ! ! output argument list: ! isli - sea/land/ice mask ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind,r_single use gridmod, only: nlat,nlon use mpimod, only: mpi_itype,mpi_comm_world,mype implicit none ! Declare passed variables integer(i_kind), intent(in ) :: iope integer(i_kind), dimension(nlat,nlon), intent( out) :: isli_anl ! Declare local variables integer(i_kind):: iret,npts !----------------------------------------------------------------------------- ! Read surface file on processor iope if(mype == iope)then call read_sfc_anl_(isli_anl) write(*,*) 'read_sfc nemsio' end if ! Load onto all processors npts=nlat*nlon call mpi_bcast(isli_anl,npts,mpi_itype,iope,mpi_comm_world,iret) end subroutine read_nemssfc_anl_ subroutine read_hsst_(hsst,cwoption) ! ! abstract : read sst bg error correlation length from berror fix file ! ! 03/13/2020 - Xu Li (based on read_wgt) ! use kinds, only : i_kind,r_single,r_kind use gridmod, only : nlat,nlon,nsig use mpeu_util, only : check_iostat implicit none integer(i_kind) , intent(in ) :: cwoption real(r_single),dimension(nlat,nlon), intent(out ) :: hsst real(r_single),dimension(nlat,nlon) :: corsst character(len=*),parameter :: myname_=myname//'::read_hsst' real(r_single),dimension(nlat,nsig,nsig):: agvin real(r_single),dimension(nlat,nsig) :: wgvin,bvin integer(i_kind) :: inerr,istat,ier integer(i_kind) :: nsigstat,nlatstat integer(i_kind) :: isig character(len=5) :: var real(r_single), allocatable, dimension(:,:) :: hwllin real(r_single), allocatable, dimension(:,:) :: corzin real(r_single), allocatable, dimension(:,:) :: corq2 real(r_single), allocatable, dimension(:,:) :: vscalesin character(len=256) :: berror_stats = "berror_stats" ! filename ! Open background error statistics file inerr = 22 open(inerr,file=berror_stats,form='unformatted',status='old',iostat=ier) call check_iostat(ier,myname_,'open('//trim(berror_stats)//')') ! Read header. Ensure that vertical resolution is consistent ! with that specified via the user namelist rewind inerr read(inerr,iostat=ier)nsigstat,nlatstat call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (nsigstat,nlatstat)') if ( nsigstat/=nsig .or. nlatstat/=nlat ) then write(6,*)'PREBAL: **ERROR** resolution of berror_stats incompatiable with GSI' write(6,*)'PREBAL: berror nsigstat,nlatstat=', nsigstat,nlatstat, & ' -vs- GSI nsig,nlat=',nsig,nlat call stop2(101) endif write(6,*) myname_,'(read_hsst): read error amplitudes ', & '"',trim(berror_stats),'". ', & 'nsigstat,nlatstat =',nsigstat,nlatstat read(inerr,iostat=ier) agvin,bvin,wgvin call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (agvin,bvin,wgvin)') write(*,*) 'in read_hsst_, cwoption = ',cwoption readloop: do read(inerr,iostat=istat) var, isig if ( istat/=0 ) exit write(*,*) 'var, isig, istat : ',var, isig, istat allocate(corzin(nlat,isig)) if ( var=='q' .or. var=='cw' ) allocate(corq2(nlat,isig)) allocate(hwllin(nlat,isig)) if ( isig>1 ) allocate(vscalesin(nlat,isig)) if ( var/='sst' ) then if ( var=='q' .or. var=='Q' .or. (var=='cw' .and. cwoption==2) ) then read(inerr,iostat=ier) corzin,corq2 call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corzin,corq2)') else read(inerr,iostat=ier) corzin call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corzin)') endif read(inerr,iostat=ier) hwllin call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (hwllin)') if ( isig>1 ) then read(inerr,iostat=ier) vscalesin call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (vscalein)') endif else read(inerr,iostat=ier) corsst call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corsst)') read(inerr,iostat=ier) hsst call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (hsst)') endif deallocate(corzin,hwllin) if (isig>1) deallocate(vscalesin) if ( var=='q' .or. var=='cw' ) deallocate(corq2) enddo readloop close(inerr) return end subroutine read_hsst_ subroutine berror_read_hsst_(iope,hsst) !$$$ subprogram documentation block ! . . . . ! subprogram: berror_read_hsst_ read hsst from berror_stats file ! prgmmr: xuli org: np23 date: 2020-03-13 ! ! program history log: ! ! input argument list: ! iope - mpi task handling i/o ! ! output argument list: ! hsst - sst bg error correlation length (km) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind,r_single use gridmod, only: nlat,nlon use mpimod, only: mpi_itype,mpi_rtype4,mpi_comm_world,mype use jfunc, only: cwoption implicit none ! Declare passed variables integer(i_kind), intent(in) :: iope real(r_single), dimension(nlat,nlon), intent(out) :: hsst ! Declare local variables integer(i_kind):: iret,npts !----------------------------------------------------------------------------- ! Read surface history file on processor iope if(mype == iope)then call read_hsst_(hsst,cwoption) write(*,*) 'read hsst from berror_stats' endif ! Load onto all processors npts=nlat*nlon call mpi_bcast(hsst,npts,mpi_rtype4,iope,mpi_comm_world,iret) end subroutine berror_read_hsst_ subroutine read_nst_ (tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d) !$$$ subprogram documentation block ! . . . . ! subprogram: read_nst_ read nems nst surface guess file (quadratic ! Gaussin grids) without scattering to tasks ! prgmmr: Huang org: np23 date: 2011-11-01 ! ! abstract: read nems surface NST file ! ! program history log: ! 2011-11-01 Huang Initial version based on sub read_gfsnst ! 2016-03-13 Li Modify for more effective I/O ! ! input argument list: ! ! output argument list: ! tref (:,:) ! oceanic foundation temperature ! dt_cool (:,:) ! sub-layer cooling amount at sub-skin layer ! z_c (:,:) ! depth of sub-layer cooling layer ! dt_warm (:,:) ! diurnal warming amount at sea surface (skin layer) ! z_w (:,:) ! depth of diurnal warming layer ! c_0 (:,:) ! coefficient to calculate d(Tz)/d(tr) ! c_d (:,:) ! coefficient to calculate d(Tz)/d(tr) ! w_0 (:,:) ! coefficient to calculate d(Tz)/d(tr) ! w_d (:,:) ! coefficient to calculate d(Tz)/d(tr) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind,r_single use mpimod, only: mype use gridmod, only: nlat_sfc,nlon_sfc use constants, only: zero,two use guess_grids, only: nfldnst,ifilenst use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close use nemsio_module, only: nemsio_gfile,nemsio_getfilehead,nemsio_readrecv implicit none ! Declare passed variables real(r_single) , dimension(nlat_sfc,nlon_sfc,nfldnst), intent( out) :: & tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d ! Declare local parameters integer(i_kind),parameter :: n_nst=9 integer(i_kind),dimension(7) :: idate integer(i_kind),dimension(4) :: odate ! Declare local variables character(len=6) :: filename character(len=120) :: my_name = 'READ_NEMSNST' character(len=1) :: null = ' ' integer(i_kind) :: i,j,it,latb,lonb integer(i_kind) :: iret, nframe integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd integer(i_kind) :: istop = 103 real(r_single) :: fhour real(r_single), dimension(nlat_sfc,nlon_sfc,nfldnst) :: xt real(r_single), allocatable, dimension(:) :: rwork2d real(r_single), allocatable, dimension(:,:) :: work !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Define read variable property !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! type(nemsio_gfile) :: gfile !----------------------------------------------------------------------------- call nemsio_init(iret=iret) if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret) do it=1,nfldnst ! read a nst file on the task write(filename,200)ifilenst(it) 200 format('nstf',i2.2) call nemsio_open(gfile,trim(filename),'READ',iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret) call nemsio_getfilehead(gfile, idate=idate, iret=iret, nframe=nframe, & nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, & dimx=lonb, dimy=latb ) if( nframe /= 0 ) then if ( mype == 0 ) & write(6,*)trim(my_name),': ***ERROR*** nframe /= 0 for global model read, nframe = ', nframe call stop2(istop) end if fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600 odate(1) = idate(4) !hour odate(2) = idate(2) !month odate(3) = idate(3) !day odate(4) = idate(1) !year if ( (latb /= nlat_sfc-2) .or. (lonb /= nlon_sfc) ) then if ( mype == 0 ) & write(6,'(a,'': inconsistent spatial dimension nlon,nlatm2 = '',2(i4,tr1),''-vs- sfc file lonb,latb = '',i4)') & trim(my_name),nlon_sfc,nlat_sfc-2,lonb,latb call stop2(80) endif ! ! Load surface fields into local work array ! Follow NEMS/GFS sfcf read order ! allocate(work(lonb,latb)) allocate(rwork2d(size(work,1)*size(work,2))) work = zero rwork2d = zero ! Tref call nemsio_readrecv(gfile, 'tref', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tref','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,tref(1,1,it),lonb,latb) ! dt_cool call nemsio_readrecv(gfile, 'dtcool','sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dt_cool','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,dt_cool(1,1,it),lonb,latb) ! z_c call nemsio_readrecv(gfile, 'zc', 'sfc', 1, rwork2d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'z_c','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,z_c(1,1,it),lonb,latb) ! xt call nemsio_readrecv(gfile, 'xt', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'xt','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,xt(1,1,it),lonb,latb) ! xz call nemsio_readrecv(gfile, 'xz', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'xz','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,z_w(1,1,it),lonb,latb) ! ! c_0 call nemsio_readrecv(gfile, 'c0', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'c_0','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,c_0(1,1,it),lonb,latb) ! c_d call nemsio_readrecv(gfile, 'cd', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'c_d','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,c_d(1,1,it),lonb,latb) ! w_0 call nemsio_readrecv(gfile, 'w0', 'sfc', 1, rwork2d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'w_0','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,w_0(1,1,it),lonb,latb) ! w_d call nemsio_readrecv(gfile, 'wd', 'sfc', 1, rwork2d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'w_d','read',istop,iret) work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/)) call tran_gfssfc(work,w_d(1,1,it),lonb,latb) ! ! Get diurnal warming amout at z=0 ! do j = 1,nlon_sfc do i = 1,nlat_sfc if (z_w(i,j,it) > zero) then dt_warm(i,j,it) = two*xt(i,j,it)/z_w(i,j,it) end if end do end do ! Deallocate local work arrays deallocate(work,rwork2d) call nemsio_close(gfile,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret) ! End of loop over time levels end do end subroutine read_nst_ subroutine read_nemsnst_ (iope,tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d) !$$$ subprogram documentation block ! . . . . ! subprogram: read_nems_nst ! prgmmr: li org: np23 date: 2016-03-13 ! ! abstract: read nems nst fields from a specific task and then broadcast to others ! ! input argument list: ! iope - mpi task handling i/o ! ! output argument list: ! tref (:,:) ! oceanic foundation temperature ! dt_cool (:,:) ! sub-layer cooling amount at sub-skin layer ! z_c (:,:) ! depth of sub-layer cooling layer ! dt_warm (:,:) ! diurnal warming amount at sea surface ! z_w (:,:) ! depth of diurnal warming layer ! c_0 (:,:) ! coefficient to calculate d(Tz)/d(tr) in dimensionless ! c_d (:,:) ! coefficient to calculate d(Tz)/d(tr) in m^-1 ! w_0 (:,:) ! coefficient to calculate d(Tz)/d(tr) in dimensionless ! w_d (:,:) ! coefficient to calculate d(Tz)/d(tr) in m^-1 ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind,r_single use gridmod, only: nlat_sfc,nlon_sfc use guess_grids, only: nfldnst use mpimod, only: mpi_itype,mpi_rtype4,mpi_comm_world use mpimod, only: mype use constants, only: zero implicit none ! Declare passed variables integer(i_kind), intent(in ) :: iope real(r_single), dimension(nlat_sfc,nlon_sfc,nfldnst), intent( out) :: & tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d ! Declare local variables integer(i_kind):: iret,npts,nptsall !----------------------------------------------------------------------------- ! Read nst file on processor iope if(mype == iope)then write(*,*) 'read_nst nemsio' call read_nst_(tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d) end if ! Load onto all processors npts=nlat_sfc*nlon_sfc nptsall=npts*nfldnst call mpi_bcast(tref, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(dt_cool, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(z_c, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(dt_warm, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(z_w, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(c_0, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(c_d, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(w_0, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) call mpi_bcast(w_d, nptsall,mpi_rtype4,iope,mpi_comm_world,iret) end subroutine read_nemsnst_ subroutine write_fv3atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) !$$$ subprogram documentation block ! . . . ! subprogram: write_fv3atm --- Gather, transform, and write out ! ! prgmmr: eliu org: np23 date: 2018-05-15 ! ! abstract: This routine gathers fields needed for the GSI analysis ! file from subdomains and then transforms the fields from ! analysis grid to model guess grid, then written to an ! atmospheric analysis file. ! ! program history log: ! 2018-05-19 eliu Initial version. Based on write_nemsatm (Huang) ! ! input argument list: ! filename - file to open and write to ! mype_out - mpi task to write output file ! gfs_bundle - bundle containing fields on subdomains ! ibin - time bin ! gfschem_bundle - (optional) bundle containing chemistry fields ! ! output argument list: ! ! attributes: ! language: f90 ! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP ! !$$$ end documentation block ! !USES: use kinds, only: r_kind,i_kind use constants, only: r1000,fv,one,zero,qcmin,r0_05,t0c use mpimod, only: mpi_rtype use mpimod, only: mpi_comm_world use mpimod, only: ierror use mpimod, only: mype use guess_grids, only: ifilesig use guess_grids, only: load_geop_hgt,geop_hgti,ges_geopi use gridmod, only: ntracer use gridmod, only: ncloud use gridmod, only: strip,jcap_b,bk5 use general_commvars_mod, only: load_grid,fill2_ns,filluv2_ns use general_specmod, only: spec_vars use obsmod, only: iadate use nemsio_module, only: nemsio_gfile,nemsio_open,nemsio_init,& nemsio_getfilehead,nemsio_close,nemsio_writerecv,nemsio_readrecv use gsi_4dvar, only: ibdate,nhr_obsbin,lwrite4danl use general_sub2grid_mod, only: sub2grid_info use egrid2agrid_mod,only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid use constants, only: two,pi,half,deg2rad use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use control_vectors, only: cvars3d use cloud_efr_mod, only: cloud_calc_gfs use mpeu_util, only: getindex implicit none ! !INPUT PARAMETERS: type(sub2grid_info), intent(in) :: grd type(spec_vars), intent(in) :: sp_a character(len=24), intent(in) :: filename ! file to open and write to integer(i_kind), intent(in) :: mype_out ! mpi task to write output file type(gsi_bundle), intent(in) :: gfs_bundle integer(i_kind), intent(in) :: ibin ! time bin !------------------------------------------------------------------------- real(r_kind),parameter:: r0_001 = 0.001_r_kind character(6):: fname_ges character(len=120) :: my_name = 'WRITE_FV3ATM_NEMS' character(len=1) :: null = ' ' integer(i_kind),dimension(7):: idate, jdate integer(i_kind),dimension(4):: odate integer(i_kind) :: k, mm1, nlatm2, nord_int, i, j, kk integer(i_kind) :: iret, lonb, latb, levs, istatus integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd integer(i_kind) :: istop = 104 integer(i_kind),dimension(5):: mydate integer(i_kind),dimension(8) :: ida,jda real(r_kind),dimension(5) :: fha real(r_kind) :: fhour real(r_kind),pointer,dimension(:,:) :: sub_ps real(r_kind),pointer,dimension(:,:,:) :: sub_u,sub_v,sub_tv real(r_kind),pointer,dimension(:,:,:) :: sub_q,sub_oz real(r_kind),pointer,dimension(:,:,:) :: sub_ql,sub_qi,sub_qr,sub_qs,sub_qg real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig) :: sub_dzb,sub_dza real(r_kind),dimension(grd%lat1*grd%lon1) :: psm real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig):: sub_dp real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: tvsm, usm, vsm real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: dpsm, qsm, ozsm real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: dzsm real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: qlsm,qism,qrsm,qssm,qgsm real(r_kind),dimension(max(grd%iglobal,grd%itotsub)) :: work1,work2 real(r_kind),dimension(grd%nlon,grd%nlat-2):: grid real(r_kind),allocatable,dimension(:) :: rwork1d,rwork1d1,rlats,rlons,clons,slons real(4),allocatable,dimension(:) :: r4lats,r4lons real(r_kind),allocatable,dimension(:,:) :: grid_b,grid_b2 real(r_kind),allocatable,dimension(:,:,:) :: grid_c, grid3, grid_c2 type(nemsio_gfile) :: gfile,gfileo logical diff_res,eqspace logical lql,lqi,lqr,lqs,lqg logical,dimension(1) :: vector type(egrid2agrid_parm) :: p_low,p_high !************************************************************************* ! Initialize local variables mm1=mype+1 nlatm2=grd%nlat-2 diff_res=.false. lql=getindex(cvars3d,'ql')>0 lqi=getindex(cvars3d,'qi')>0 lqr=getindex(cvars3d,'qr')>0 lqs=getindex(cvars3d,'qs')>0 lqg=getindex(cvars3d,'qg')>0 istatus=0 call gsi_bundlegetpointer(gfs_bundle,'ps', sub_ps,iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'u', sub_u, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'v', sub_v, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'tv', sub_tv,iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'q', sub_q, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'oz', sub_oz,iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'ql', sub_ql,iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'qi', sub_qi,iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'qr', sub_qr,iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'qs', sub_qs,iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'qg', sub_qg,iret); istatus=istatus+iret if ( istatus /= 0 ) then if ( mype == 0 ) then write(6,*) 'write_fv3atm_: ERROR' write(6,*) 'Missing some of the required fields' write(6,*) 'Aborting ... ' endif call stop2(999) endif if ( sp_a%jcap /= jcap_b ) then if ( mype == 0 ) write(6, & '('' dual resolution for nems sp_a%jcap,jcap_b = '',2i6)') & sp_a%jcap,jcap_b diff_res = .true. endif ! Single task writes analysis data to analysis file if ( mype == mype_out ) then write(fname_ges,'(''sigf'',i2.2)') ifilesig(ibin) ! Read header information from first guess file. call nemsio_init(iret) if ( iret /= 0 ) call error_msg(trim(my_name),null,null,'init',istop,iret) call nemsio_open(gfile,trim(fname_ges),'read',iret) if ( iret /= 0 ) call error_msg(trim(my_name),trim(fname_ges),null,'open',istop,iret) call nemsio_getfilehead(gfile, iret=iret, nfhour=nfhour, & nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, & idate=idate, dimx=lonb, dimy=latb, dimz=levs) if ( iret /= 0 ) then write(6,*) trim(my_name),': problem with nemsio_getfilehead, Status = ',iret call stop2(103) endif if ( levs /= grd%nsig ) then write(6,*) trim(my_name),': problem in data dimension background levs = ',levs,' nsig = ',grd%nsig call stop2(103) endif ! copy input header info to output header info gfileo=gfile ! Update header information (with ibdate) and write it to analysis file ! (w/ _open statement). mydate=ibdate fha(:)=zero ; ida=0; jda=0 fha(2)=real(nhr_obsbin*(ibin-1)) ! relative time interval in hours ida(1)=mydate(1) ! year ida(2)=mydate(2) ! month ida(3)=mydate(3) ! day ida(4)=0 ! time zone ida(5)=mydate(4) ! hour ! Move date-time forward by nhr_assimilation hours call w3movdat(fha,ida,jda) jdate(1) = jda(1) ! analysis year jdate(2) = jda(2) ! analysis month jdate(3) = jda(3) ! analysis day jdate(4) = jda(5) ! analysis hour jdate(5) = iadate(5) ! analysis minute jdate(6) = 0 ! analysis scaled seconds jdate(7) = idate(7) ! analysis seconds multiplier nfhour =0 ! new forecast hour, zero at analysis time nfminute =0 nfsecondn=0 nfsecondd=100 ! default for denominator fhour = zero odate(1) = jdate(4) !hour odate(2) = jdate(2) !month odate(3) = jdate(3) !day odate(4) = jdate(1) !year ! open new output file with new header gfileo with "write" access. ! Use this call to update header as well call nemsio_open(gfileo,trim(filename),'write',iret=iret, & idate=jdate, nfhour=nfhour, nfminute=nfminute, & nfsecondn=nfsecondn, nfsecondd=nfsecondd) if ( iret /= 0 ) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret) ! Allocate structure arrays to hold data allocate(rwork1d(latb*lonb),rwork1d1(latb*lonb)) if ( diff_res .or. lupp) then allocate( grid_b(lonb,latb),grid_c(latb+2,lonb,1),grid3(grd%nlat,grd%nlon,1)) allocate( grid_b2(lonb,latb),grid_c2(latb+2,lonb,1)) allocate( rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb),r4lats(lonb*latb),r4lons(lonb*latb)) call nemsio_getfilehead(gfile,lat=r4lats,iret=iret) call nemsio_getfilehead(gfile,lon=r4lons,iret=iret) do j=1,latb rlats(latb+2-j)=deg2rad*r4lats(lonb/2+(j-1)*lonb) enddo rlats(1)=-half*pi rlats(latb+2)=half*pi do j=1,lonb rlons(j)=deg2rad*r4lons(j) enddo do j=1,lonb clons(j)=cos(rlons(j)) slons(j)=sin(rlons(j)) enddo nord_int=4 eqspace=.false. call g_create_egrid2agrid(grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons, & latb+2,rlats,lonb,rlons,& nord_int,p_low,.false.,eqspace=eqspace) call g_create_egrid2agrid(latb+2,rlats,lonb,rlons, & grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons,& nord_int,p_high,.false.,eqspace=eqspace) deallocate(rlats,rlons,r4lats,r4lons) endif ! if ( diff_res ) ! Terrain ! Write out input file surface height call nemsio_readrecv(gfile,'hgt', 'sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'hgt','writeread',istop,iret) call nemsio_writerecv(gfileo,'hgt','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'hgt','write',istop,iret) endif ! if ( mype == mype_out ) ! Calculate delz increment for UPP if (lupp) then do k=1,grd%nsig ! k goes from bottom to top, so k - k+1 is negative delz sub_dzb(:,:,k) = ges_geopi(:,:,k,ibin) - ges_geopi(:,:,k+1,ibin) enddo if ((.not. lwrite4danl) .or. ibin == 1) call load_geop_hgt do k=1,grd%nsig sub_dza(:,:,k) = geop_hgti(:,:,k,ibin) - geop_hgti(:,:,k+1,ibin) enddo sub_dza = sub_dza - sub_dzb !sub_dza is increment endif ! Strip off boundary points from subdomains call strip(sub_ps ,psm) call strip(sub_tv ,tvsm ,grd%nsig) call strip(sub_q ,qsm ,grd%nsig) call strip(sub_oz ,ozsm ,grd%nsig) call strip(sub_dp ,dpsm ,grd%nsig) call strip(sub_u ,usm ,grd%nsig) call strip(sub_v ,vsm ,grd%nsig) if (lql ) call strip(sub_ql ,qlsm ,grd%nsig) if (lqi ) call strip(sub_qi ,qism ,grd%nsig) if (lqr ) call strip(sub_qr ,qrsm ,grd%nsig) if (lqs ) call strip(sub_qs ,qssm ,grd%nsig) if (lqg ) call strip(sub_qg ,qgsm ,grd%nsig) if (lupp) call strip(sub_dza ,dzsm ,grd%nsig) ! Thermodynamic variable ! The GSI analysis variable is virtual temperature (Tv). For NEMSIO ! output we need the sensible temperature. ! Convert Tv to T tvsm = tvsm/(one+fv*qsm) ! Generate and write analysis fields ! Surface pressure. call mpi_gatherv(psm,grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype==mype_out) then if(diff_res .or. lupp)then call nemsio_readrecv(gfile,'pres','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pres','read',istop,iret) rwork1d1 = r0_001*rwork1d grid_b=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) if (lupp) work1(kk)=grid3(i,j,1) end do if (lupp) then do k=1,grd%nsig do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)*(bk5(k)-bk5(k+1)) enddo call g_egrid2agrid(p_high,grid3,grid_c2,1,1,vector) call nemsio_readrecv(gfile,'dpres','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dpres','read',istop,iret) grid_b2=reshape(rwork1d,(/size(grid_b2,1),size(grid_b2,2)/)) do j=1,latb do i=1,lonb grid_b2(i,j)=grid_b2(i,j)+r1000*(grid_c2(latb-j+2,i,1)) enddo enddo rwork1d = reshape(grid_b2,(/size(rwork1d)/)) call nemsio_writerecv(gfileo,'dpres','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dpres','write',istop,iret) enddo do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk) enddo endif call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=r1000*(grid_b(i,j)+grid_c(latb-j+2,i,1)) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) grid = grid*r1000 rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'pres','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'psfc','write',istop,iret) endif ! u, v do k=1,grd%nsig call mpi_gatherv(usm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) call mpi_gatherv(vsm(1,k),grd%ijn(mm1),mpi_rtype,& work2,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype==mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'ugrd','mid layer',k,rwork1d,iret=iret) call nemsio_readrecv(gfile,'vgrd','mid layer',k,rwork1d1,iret=iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) grid_b2=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.true. call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do call g_egrid2agrid(p_low,grid_c2,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work2(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b2(i,j)=grid_b2(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) rwork1d1 = reshape(grid_b2,(/size(rwork1d1)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) call load_grid(work2,grid) rwork1d1 = reshape(grid,(/size(rwork1d1)/)) end if ! Zonal wind call nemsio_writerecv(gfileo,'ugrd','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ugrd','write',istop,iret) ! Meridional wind call nemsio_writerecv(gfileo,'vgrd','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vgrd','write',istop,iret) endif end do ! Thermodynamic variable do k=1,grd%nsig call mpi_gatherv(tvsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'tmp','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'tmp','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','write',istop,iret) endif end do ! Specific humidity do k=1,grd%nsig call mpi_gatherv(qsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'spfh','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'spfh','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','write',istop,iret) endif end do ! Ozone do k=1,grd%nsig call mpi_gatherv(ozsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'o3mr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'o3mr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','write',istop,iret) endif end do ! Cloud condensate mixing ratio if (ntracer>2 .or. ncloud>=1) then ! Cloud liquid water do k=1,grd%nsig if (lql) then call mpi_gatherv(qlsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'clwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-max(grid3(i,j,1),qcmin) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'clwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','write',istop,iret) endif else if (mype == mype_out) then call nemsio_readrecv(gfile,'clwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','read',istop,iret) call nemsio_writerecv(gfileo,'clwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','write',istop,iret) endif endif end do ! Cloud ice do k=1,grd%nsig if (lqi) then call mpi_gatherv(qism(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'icmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-max(grid3(i,j,1),qcmin) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'icmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','write',istop,iret) endif else if (mype == mype_out) then call nemsio_readrecv(gfile,'icmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','read',istop,iret) call nemsio_writerecv(gfileo,'icmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','write',istop,iret) endif endif end do ! Rain do k=1,grd%nsig if (lqr) then call mpi_gatherv(qrsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'rwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'rwmr','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-max(grid3(i,j,1),qcmin) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'rwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'rwmr','write',istop,iret) endif else if (mype == mype_out) then call nemsio_readrecv(gfile,'rwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'rwmr','read',istop,iret) call nemsio_writerecv(gfileo,'rwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'rwmr','write',istop,iret) endif endif end do ! Snow do k=1,grd%nsig if (lqs) then call mpi_gatherv(qssm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'snmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'snmr','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-max(grid3(i,j,1),qcmin) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'snmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'snmr','write',istop,iret) endif else if (mype == mype_out) then call nemsio_readrecv(gfile,'snmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'snmr','read',istop,iret) call nemsio_writerecv(gfileo,'snmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'snmr','write',istop,iret) endif endif end do ! Graupel do k=1,grd%nsig if (lqg) then call mpi_gatherv(qgsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'grle','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'grle','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-max(grid3(i,j,1),qcmin) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'grle','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'grle','write',istop,iret) endif else if (mype == mype_out) then call nemsio_readrecv(gfile,'grle','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'grle','read',istop,iret) call nemsio_writerecv(gfileo,'grle','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'grle','write',istop,iret) endif endif end do ! Cloud Amount (cloud fraction) - should be the same as the input guess values do k=1,grd%nsig if (mype==mype_out) then call nemsio_readrecv(gfile,'cld_amt','mid layer',k,rwork1d,iret=iret) if (iret == 0) then call nemsio_writerecv(gfileo,'cld_amt','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'cld_amt','write',istop,iret) endif endif enddo endif !ntracer ! Variables needed by the Unified Post Processor (dzdt, delz, delp) if (lupp) then if (mype == mype_out) then do k=1,grd%nsig call nemsio_readrecv(gfile,'dzdt','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dzdt','read',istop,iret) call nemsio_writerecv(gfileo,'dzdt','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dzdt','write',istop,iret) enddo endif do k=1,grd%nsig call mpi_gatherv(dzsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then call nemsio_readrecv(gfile,'delz','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','read',istop,iret) ! if delz in background is positive, flip sign of increment if (sum(rwork1d) > zero) work1 = -one * work1 if(diff_res)then grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = rwork1d + reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'delz','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','write',istop,iret) endif end do endif ! ! Deallocate local array ! if (mype==mype_out) then if (diff_res .or. lupp) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid3,clons,slons) call nemsio_close(gfile,iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),null,'close',istop,iret) call nemsio_close(gfileo,iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret) ! ! Deallocate local array ! deallocate(rwork1d,rwork1d1) ! write(6,'(a,'': atm anal written for lonb,latb,levs= '',3i6,'',valid hour= '',f4.1,'',idate= '',4i5)') & trim(my_name),lonb,latb,levs,fhour,odate endif end subroutine write_fv3atm_ subroutine write_atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin,gfschem_bundle) !$$$ subprogram documentation block ! . . . ! subprogram: write_nemsatm --- Gather, transform, and write out ! ! prgmmr: Huang org: np23 date: 2010-02-22 ! ! abstract: This routine gathers fields needed for the GSI analysis ! file from subdomains and then transforms the fields from ! analysis grid to model guess grid, then written to an ! atmospheric analysis file. ! ! program history log: ! 2010-02-22 Huang Initial version. Based on write_gfsatm ! 2011-02-14 Huang Re-arrange the write sequence to be same as model ! read/rite sequence. ! 2013-10-25 todling reposition load_grid to commvars ! 2016-07-28 mahajan update with bundling ability ! 2019-04-19 Wei/Martin - added gfschem_bundle to write optional aerosols ! for both FV3-Chem and NGAC ! ! input argument list: ! filename - file to open and write to ! mype_out - mpi task to write output file ! gfs_bundle - bundle containing fields on subdomains ! ibin - time bin ! ! output argument list: ! ! attributes: ! language: f90 ! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP ! !$$$ end documentation block ! !USES: use kinds, only: r_kind,i_kind use constants, only: r1000,fv,one,zero,qcmin,r0_05,t0c use mpimod, only: mpi_rtype use mpimod, only: mpi_comm_world use mpimod, only: ierror use mpimod, only: mype use guess_grids, only: ifilesig use guess_grids, only: ges_prsi use guess_grids, only: load_geop_hgt,geop_hgti,ges_geopi use gridmod, only: ntracer use gridmod, only: ncloud use gridmod, only: strip,jcap_b,bk5 use gridmod, only: use_fv3_aero use general_commvars_mod, only: load_grid,fill2_ns,filluv2_ns use general_specmod, only: spec_vars use obsmod, only: iadate use nemsio_module, only: nemsio_gfile,nemsio_open,nemsio_init,& nemsio_getfilehead,nemsio_close,nemsio_writerecv,nemsio_readrecv use gsi_4dvar, only: ibdate,nhr_obsbin,lwrite4danl use general_sub2grid_mod, only: sub2grid_info use egrid2agrid_mod,only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid use constants, only: two,pi,half,deg2rad use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use cloud_efr_mod, only: cloud_calc_gfs use chemmod, only: laeroana_gocart implicit none ! !INPUT PARAMETERS: type(sub2grid_info), intent(in) :: grd type(spec_vars), intent(in) :: sp_a character(len=24), intent(in) :: filename ! file to open and write to integer(i_kind), intent(in) :: mype_out ! mpi task to write output file type(gsi_bundle), intent(in) :: gfs_bundle type(gsi_bundle),optional,intent(in) :: gfschem_bundle ! for aerosols integer(i_kind), intent(in) :: ibin ! time bin !------------------------------------------------------------------------- real(r_kind),parameter:: r0_001 = 0.001_r_kind character(6):: fname_ges character(len=120) :: my_name = 'WRITE_NEMSATM' character(len=1) :: null = ' ' integer(i_kind),dimension(7):: idate, jdate integer(i_kind),dimension(4):: odate integer(i_kind) :: k, mm1, nlatm2, nord_int, i, j, kk integer(i_kind) :: iret, lonb, latb, levs, istatus integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd integer(i_kind) :: istop = 104 integer(i_kind),dimension(5):: mydate integer(i_kind),dimension(8) :: ida,jda real(r_kind),dimension(5) :: fha real(r_kind) :: fhour real(r_kind),pointer,dimension(:,:) :: sub_ps real(r_kind),pointer,dimension(:,:,:) :: sub_u,sub_v,sub_tv real(r_kind),pointer,dimension(:,:,:) :: sub_q,sub_oz,sub_cwmr ! Sub-domain aerosol arrays real(r_kind),pointer,dimension(:,:,:) :: sub_du1,sub_du2,sub_du3,sub_du4,sub_du5 real(r_kind),pointer,dimension(:,:,:) :: sub_ss1,sub_ss2,sub_ss3,sub_ss4,sub_so4 real(r_kind),pointer,dimension(:,:,:) :: sub_oc1,sub_oc2,sub_bc1,sub_bc2 real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig) :: sub_dzb,sub_dza real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig+1) :: sub_prsi real(r_kind),dimension(grd%lat1*grd%lon1) :: psm real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig):: sub_dp real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: tvsm, usm, vsm real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: dpsm, qsm, ozsm real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: cwsm, dzsm ! Aerosol array real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: & du001sm,du002sm,du003sm,du004sm,du005sm, & ss001sm,ss002sm,ss003sm,ss004sm, & so4sm,ocphosm,ocphism,bcphosm,bcphism integer(i_kind) :: m real(r_kind),dimension(max(grd%iglobal,grd%itotsub)) :: work1,work2 real(r_kind),dimension(grd%nlon,grd%nlat-2):: grid real(r_kind),allocatable,dimension(:) :: rwork1d,rwork1d1,rlats,rlons,clons,slons real(4),allocatable,dimension(:) :: r4lats,r4lons real(r_kind),allocatable,dimension(:,:) :: grid_b,grid_b2 real(r_kind),allocatable,dimension(:,:,:) :: grid_c, grid3, grid_c2, grid3b type(nemsio_gfile) :: gfile,gfileo logical diff_res,eqspace logical,dimension(1) :: vector type(egrid2agrid_parm) :: p_low,p_high !************************************************************************* ! Initialize local variables mm1=mype+1 nlatm2=grd%nlat-2 diff_res=.false. istatus=0 call gsi_bundlegetpointer(gfs_bundle,'ps', sub_ps, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'u', sub_u, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'v', sub_v, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'tv', sub_tv, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'q', sub_q, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'oz', sub_oz, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfs_bundle,'cw', sub_cwmr,iret); istatus=istatus+iret if ( istatus /= 0 ) then if ( mype == 0 ) then write(6,*) 'write_atm_: ERROR' write(6,*) 'Missing some of the required fields' write(6,*) 'Aborting ... ' endif call stop2(999) endif if (present(gfschem_bundle) .and. laeroana_gocart) then call gsi_bundlegetpointer(gfschem_bundle,'sulf', sub_so4, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'oc1', sub_oc1, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'oc2', sub_oc2, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'bc1', sub_bc1, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'bc2', sub_bc2, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'dust1', sub_du1, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'dust2', sub_du2, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'dust3', sub_du3, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'dust4', sub_du4, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'dust5', sub_du5, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'seas1', sub_ss1, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'seas2', sub_ss2, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'seas3', sub_ss3, iret); istatus=istatus+iret call gsi_bundlegetpointer(gfschem_bundle,'seas4', sub_ss4, iret); istatus=istatus+iret end if if ( sp_a%jcap /= jcap_b ) then if ( mype == 0 ) write(6, & '('' dual resolution for nems sp_a%jcap,jcap_b = '',2i6)') & sp_a%jcap,jcap_b diff_res = .true. endif ! Single task writes analysis data to analysis file if ( mype == mype_out ) then write(fname_ges,'(''sigf'',i2.2)') ifilesig(ibin) ! Read header information from first guess file. call nemsio_init(iret) if ( iret /= 0 ) call error_msg(trim(my_name),null,null,'init',istop,iret) call nemsio_open(gfile,trim(fname_ges),'read',iret) if ( iret /= 0 ) call error_msg(trim(my_name),trim(fname_ges),null,'open',istop,iret) call nemsio_getfilehead(gfile, iret=iret, nfhour=nfhour, & nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, & idate=idate, dimx=lonb, dimy=latb, dimz=levs) if ( iret /= 0 ) then write(6,*) trim(my_name),': problem with nemsio_getfilehead, Status = ',iret call stop2(103) endif if ( levs /= grd%nsig ) then write(6,*) trim(my_name),': problem in data dimension background levs = ',levs,' nsig = ',grd%nsig call stop2(103) endif ! copy input header info to output header info gfileo=gfile ! Update header information (with ibdate) and write it to analysis file (w/ _open statement). mydate=ibdate fha(:)=zero ; ida=0; jda=0 fha(2)=real(nhr_obsbin*(ibin-1)) ! relative time interval in hours ida(1)=mydate(1) ! year ida(2)=mydate(2) ! month ida(3)=mydate(3) ! day ida(4)=0 ! time zone ida(5)=mydate(4) ! hour ! Move date-time forward by nhr_assimilation hours call w3movdat(fha,ida,jda) jdate(1) = jda(1) ! analysis year jdate(2) = jda(2) ! analysis month jdate(3) = jda(3) ! analysis day jdate(4) = jda(5) ! analysis hour jdate(5) = iadate(5) ! analysis minute jdate(6) = 0 ! analysis scaled seconds jdate(7) = idate(7) ! analysis seconds multiplier nfhour =0 ! new forecast hour, zero at analysis time nfminute =0 nfsecondn=0 nfsecondd=100 ! default for denominator fhour = zero odate(1) = jdate(4) !hour odate(2) = jdate(2) !month odate(3) = jdate(3) !day odate(4) = jdate(1) !year ! open new output file with new header gfileo with "write" access. ! Use this call to update header as well call nemsio_open(gfileo,trim(filename),'write',iret=iret, & idate=jdate, nfhour=nfhour, nfminute=nfminute, & nfsecondn=nfsecondn, nfsecondd=nfsecondd) if ( iret /= 0 ) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret) ! Allocate structure arrays to hold data allocate(rwork1d(latb*lonb),rwork1d1(latb*lonb)) if (imp_physics == 11) allocate(grid3b(grd%nlat,grd%nlon,1)) if ( diff_res .or. imp_physics == 11 .or. lupp) then allocate( grid_b(lonb,latb),grid_c(latb+2,lonb,1),grid3(grd%nlat,grd%nlon,1)) allocate( grid_b2(lonb,latb),grid_c2(latb+2,lonb,1)) allocate( rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb),r4lats(lonb*latb),r4lons(lonb*latb)) call nemsio_getfilehead(gfile,lat=r4lats,iret=iret) call nemsio_getfilehead(gfile,lon=r4lons,iret=iret) do j=1,latb rlats(latb+2-j)=deg2rad*r4lats(lonb/2+(j-1)*lonb) enddo rlats(1)=-half*pi rlats(latb+2)=half*pi do j=1,lonb rlons(j)=deg2rad*r4lons(j) enddo do j=1,lonb clons(j)=cos(rlons(j)) slons(j)=sin(rlons(j)) enddo nord_int=4 eqspace=.false. call g_create_egrid2agrid(grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons, & latb+2,rlats,lonb,rlons,& nord_int,p_low,.false.,eqspace=eqspace) call g_create_egrid2agrid(latb+2,rlats,lonb,rlons, & grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons,& nord_int,p_high,.false.,eqspace=eqspace) deallocate(rlats,rlons,r4lats,r4lons) endif ! if ( diff_res ) ! Terrain ! Write out input file surface height call nemsio_readrecv(gfile,'hgt', 'sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'hgt','writeread',istop,iret) call nemsio_writerecv(gfileo,'hgt','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'hgt','write',istop,iret) endif ! if ( mype == mype_out ) sub_prsi = ges_prsi(:,:,:,ibin) do k=1,grd%nsig sub_dp(:,:,k) = sub_prsi(:,:,k) - sub_prsi(:,:,k+1) end do ! Calculate delz increment for UPP if (lupp) then do k=1,grd%nsig ! k goes from bottom to top, so k - k+1 is negative delz sub_dzb(:,:,k) = ges_geopi(:,:,k,ibin) - ges_geopi(:,:,k+1,ibin) enddo if ((.not. lwrite4danl) .or. ibin == 1) call load_geop_hgt do k=1,grd%nsig sub_dza(:,:,k) = geop_hgti(:,:,k,ibin) - geop_hgti(:,:,k+1,ibin) enddo sub_dza = sub_dza - sub_dzb !sub_dza is increment endif ! Strip off boundary points from subdomains call strip(sub_ps ,psm) call strip(sub_tv ,tvsm ,grd%nsig) call strip(sub_q ,qsm ,grd%nsig) call strip(sub_oz ,ozsm ,grd%nsig) call strip(sub_cwmr,cwsm ,grd%nsig) call strip(sub_dp ,dpsm ,grd%nsig) call strip(sub_u ,usm ,grd%nsig) call strip(sub_v ,vsm ,grd%nsig) if (lupp) call strip(sub_dza ,dzsm ,grd%nsig) if (laeroana_gocart) then call strip(sub_du1 ,du001sm ,grd%nsig) call strip(sub_du2 ,du002sm ,grd%nsig) call strip(sub_du3 ,du003sm ,grd%nsig) call strip(sub_du4 ,du004sm ,grd%nsig) call strip(sub_du5 ,du005sm ,grd%nsig) call strip(sub_ss1 ,ss001sm ,grd%nsig) call strip(sub_ss2 ,ss002sm ,grd%nsig) call strip(sub_ss3 ,ss003sm ,grd%nsig) call strip(sub_ss4 ,ss004sm ,grd%nsig) call strip(sub_so4 ,so4sm ,grd%nsig) call strip(sub_oc1 ,ocphosm ,grd%nsig) call strip(sub_oc2 ,ocphism ,grd%nsig) call strip(sub_bc1 ,bcphosm ,grd%nsig) call strip(sub_bc2 ,bcphism ,grd%nsig) end if ! Thermodynamic variable ! The GSI analysis variable is virtual temperature (Tv). For NEMSIO ! output we need the sensible temperature. ! Convert Tv to T tvsm = tvsm/(one+fv*qsm) ! Generate and write analysis fields ! Surface pressure. call mpi_gatherv(psm,grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype==mype_out) then if(diff_res .or. lupp)then call nemsio_readrecv(gfile,'pres','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pres','read',istop,iret) rwork1d1 = r0_001*rwork1d grid_b=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) if (lupp) work1(kk)=grid3(i,j,1) end do if (lupp) then do k=1,grd%nsig do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)*(bk5(k)-bk5(k+1)) enddo call g_egrid2agrid(p_high,grid3,grid_c2,1,1,vector) call nemsio_readrecv(gfile,'dpres','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dpres','read',istop,iret) grid_b2=reshape(rwork1d,(/size(grid_b2,1),size(grid_b2,2)/)) do j=1,latb do i=1,lonb grid_b2(i,j)=grid_b2(i,j)+r1000*(grid_c2(latb-j+2,i,1)) enddo enddo rwork1d = reshape(grid_b2,(/size(rwork1d)/)) call nemsio_writerecv(gfileo,'dpres','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dpres','write',istop,iret) enddo do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk) enddo endif call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=r1000*(grid_b(i,j)+grid_c(latb-j+2,i,1)) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) grid = grid*r1000 rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'pres','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'psfc','write',istop,iret) endif ! u, v do k=1,grd%nsig call mpi_gatherv(usm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) call mpi_gatherv(vsm(1,k),grd%ijn(mm1),mpi_rtype,& work2,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype==mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'ugrd','mid layer',k,rwork1d,iret=iret) call nemsio_readrecv(gfile,'vgrd','mid layer',k,rwork1d1,iret=iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) grid_b2=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.true. call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do call g_egrid2agrid(p_low,grid_c2,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work2(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b2(i,j)=grid_b2(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) rwork1d1 = reshape(grid_b2,(/size(rwork1d1)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) call load_grid(work2,grid) rwork1d1 = reshape(grid,(/size(rwork1d1)/)) end if ! Zonal wind call nemsio_writerecv(gfileo,'ugrd','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ugrd','write',istop,iret) ! Meridional wind call nemsio_writerecv(gfileo,'vgrd','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vgrd','write',istop,iret) endif end do ! Thermodynamic variable do k=1,grd%nsig call mpi_gatherv(tvsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'tmp','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'tmp','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','write',istop,iret) endif end do ! Specific humidity do k=1,grd%nsig call mpi_gatherv(qsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'spfh','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'spfh','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','write',istop,iret) endif end do ! Ozone do k=1,grd%nsig call mpi_gatherv(ozsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then call nemsio_readrecv(gfile,'o3mr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'o3mr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','write',istop,iret) endif end do ! Cloud condensate mixing ratio if (ntracer>2 .or. ncloud>=1) then do k=1,grd%nsig call mpi_gatherv(cwsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (imp_physics == 11) then call mpi_gatherv(tvsm(1,k),grd%ijn(mm1),mpi_rtype,& work2,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) endif if (mype == mype_out) then if(diff_res .or. imp_physics == 11)then call nemsio_readrecv(gfile,'clwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) if (imp_physics == 11) then call nemsio_readrecv(gfile,'icmr','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','read',istop,iret) grid_b2=reshape(rwork1d1,(/size(grid_b2,1),size(grid_b2,2)/)) grid_b = grid_b + grid_b2 endif vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-max(grid3(i,j,1),qcmin) if (imp_physics == 11) then work2(kk) = -r0_05*(work2(kk) - t0c) work2(kk) = max(zero,work2(kk)) work2(kk) = min(one,work2(kk)) grid3b(i,j,1)=grid3(i,j,1) grid3(i,j,1)=grid3b(i,j,1)*(one - work2(kk)) endif end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) if (imp_physics == 11) grid_b = grid_b - grid_b2 do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) if (imp_physics == 11) then do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=grid3b(i,j,1)*work2(kk) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b2(i,j)=grid_b2(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d1 = reshape(grid_b2,(/size(rwork1d1)/)) endif else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) endif call nemsio_writerecv(gfileo,'clwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','write',istop,iret) if (imp_physics == 11) then call nemsio_writerecv(gfileo,'icmr','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','write',istop,iret) if (lupp) then call nemsio_readrecv(gfile,'rwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'rwmr','read',istop,iret) call nemsio_writerecv(gfileo,'rwmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'rwmr','write',istop,iret) call nemsio_readrecv(gfile,'snmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'snmr','read',istop,iret) call nemsio_writerecv(gfileo,'snmr','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'snmr','write',istop,iret) call nemsio_readrecv(gfile,'grle','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'grle','read',istop,iret) call nemsio_writerecv(gfileo,'grle','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'grle','write',istop,iret) endif call nemsio_readrecv(gfile,'cld_amt','mid layer',k,rwork1d,iret=iret) if (iret == 0) then call nemsio_writerecv(gfileo,'cld_amt','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'cld_amt','write',istop,iret) endif endif endif !mype == mype_out end do endif !ntracer ! Variables needed by the Unified Post Processor (dzdt, delz, delp) if (lupp) then if (mype == mype_out) then do k=1,grd%nsig call nemsio_readrecv(gfile,'dzdt','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dzdt','read',istop,iret) call nemsio_writerecv(gfileo,'dzdt','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dzdt','write',istop,iret) enddo endif do k=1,grd%nsig call mpi_gatherv(dzsm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then call nemsio_readrecv(gfile,'delz','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','read',istop,iret) ! if delz in background is positive, flip sign of increment if (sum(rwork1d) > zero) work1 = -one * work1 ! Flip sign, FV3 is top to bottom if(diff_res)then grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = rwork1d + reshape(grid,(/size(rwork1d)/)) end if call nemsio_writerecv(gfileo,'delz','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','write',istop,iret) endif end do endif ! aerosol output if laeroana_gocart is T ! du001 if (laeroana_gocart) then do k=1,grd%nsig call mpi_gatherv(du001sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'dust1','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'du001','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du001','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'dust1','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'du001','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du001','write',istop,iret) endif end do ! du002 do k=1,grd%nsig call mpi_gatherv(du002sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'dust2','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'du002','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du002','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'dust2','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'du002','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du002','write',istop,iret) endif end do ! du003 do k=1,grd%nsig call mpi_gatherv(du003sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'dust3','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'du003','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du003','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'dust3','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'du003','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du003','write',istop,iret) endif end do ! du004 do k=1,grd%nsig call mpi_gatherv(du004sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'dust4','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'du004','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du004','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'dust4','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'du004','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du004','write',istop,iret) endif end do ! du005 do k=1,grd%nsig call mpi_gatherv(du005sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'dust5','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'du005','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du005','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'dust5','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'du005','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'du005','write',istop,iret) endif end do ! ss001 and ss002 NOTE: It depends the ratio of ss1 and ss2 in guess to ! distribute the analysis mixing ratio do k=1,grd%nsig call mpi_gatherv(ss001sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if (use_fv3_aero) then call nemsio_readrecv(gfile,'seas1','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss001','read',istop,iret) call nemsio_readrecv(gfile,'seas2','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss002','read',istop,iret) else call nemsio_readrecv(gfile,'ss001','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss001','read',istop,iret) call nemsio_readrecv(gfile,'ss002','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss002','read',istop,iret) end if rwork1d=rwork1d+rwork1d1 ! first guess ss1+ss2, rwork1d become the original total mixing ratio do m=1,size(rwork1d) if (rwork1d(m) == zero) then ! if original total mixing ratio is zero, which mean ss1=ss2=0. Set the rwork1d1:ss2ratio to zero rwork1d1(m)=zero else rwork1d1(m)=rwork1d1(m)/rwork1d(m) ! if original total mixing ratio isn't zero, which mean ss1!=0 or ss2!=0. ! Set ! rwork1d1:ss2ratio=rwork1d1:ss2/rwork1d:total end if end do if(diff_res)then grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) ! analysis total mixing ratio else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) ! analysis total mixing ratio end if ! if there is increment in analysis but nothing originally, half-half to ! analysis sea salt 1 and 2 respectively. ! do m=1,size(rwork1d) if (rwork1d(m) /= zero .and. rwork1d1(m) == zero ) then ! if ana!=0, ratio=0. Then ss1=ss2=half total mixing ratio rwork1d1(m)=half*rwork1d(m) ! sea salt 2 else if ( rwork1d1(m) /= zero .and. rwork1d(m) == zero ) then ! if ratio!=0, ana=0. Then ss1=ss2=0. Set ss2=0. rwork1d1(m)=zero else ! else mean ana=0 ratio=0 or ana!=0 ratio!=0. ! Then ss2=ana*ratio rwork1d1(m)=rwork1d(m)*rwork1d1(m) ! sea salt 2 end if rwork1d(m)=rwork1d(m)-rwork1d1(m) ! sea salt 1 enddo if (use_fv3_aero) then call nemsio_writerecv(gfileo,'seas1','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss001','write',istop,iret) call nemsio_writerecv(gfileo,'seas2','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss002','write',istop,iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg rwork1d1=rwork1d1*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'ss001','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss001','write',istop,iret) call nemsio_writerecv(gfileo,'ss002','mid layer',k,rwork1d1,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss002','write',istop,iret) end if endif end do ! ss003 do k=1,grd%nsig call mpi_gatherv(ss002sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'seas3','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'ss003','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss003','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'seas3','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'ss003','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss003','write',istop,iret) endif end do ! ss004 do k=1,grd%nsig call mpi_gatherv(ss003sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'seas4','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'ss004','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss004','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'seas4','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'ss004','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss004','write',istop,iret) endif end do ! ss005 do k=1,grd%nsig call mpi_gatherv(ss004sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'seas5','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'ss005','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss005','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'seas5','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'ss005','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ss005','write',istop,iret) endif end do ! dms, msa, so2 do k=1,grd%nsig if (mype == mype_out) then if (use_fv3_aero) then call nemsio_readrecv(gfile,'DMS','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'dms','mid layer',k,rwork1d,iret=iret) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'DMS','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'dms','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dms','write',istop,iret) endif end do do k=1,grd%nsig if (mype == mype_out) then call nemsio_readrecv(gfile,'msa','mid layer',k,rwork1d,iret=iret) call nemsio_writerecv(gfileo,'msa','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'msa','write',istop,iret) endif end do do k=1,grd%nsig if (mype == mype_out) then call nemsio_readrecv(gfile,'so2','mid layer',k,rwork1d,iret=iret) call nemsio_writerecv(gfileo,'so2','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'so2','write',istop,iret) endif end do ! so4 do k=1,grd%nsig call mpi_gatherv(so4sm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'sulf','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'so4','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'so4','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'sulf','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'so4','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'so4','write',istop,iret) endif end do ! oc1 do k=1,grd%nsig call mpi_gatherv(ocphosm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'oc1','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'ocphobic','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ocphobic','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'oc1','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'ocphobic','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ocphobic','write',istop,iret) endif end do ! oc2 do k=1,grd%nsig call mpi_gatherv(ocphism(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'oc2','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'ocphilic','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ocphilic','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'oc2','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'ocphilic','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ocphilic','write',istop,iret) endif end do ! bc1 do k=1,grd%nsig call mpi_gatherv(bcphosm(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'bc1','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'bcphobic','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'bcphobic','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'bc1','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'bcphobic','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'bcphobic','write',istop,iret) endif end do ! bc2 do k=1,grd%nsig call mpi_gatherv(bcphism(1,k),grd%ijn(mm1),mpi_rtype,& work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then if(diff_res)then if (use_fv3_aero) then call nemsio_readrecv(gfile,'bc2','mid layer',k,rwork1d,iret=iret) else call nemsio_readrecv(gfile,'bcphilic','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'bcphilic','read',istop,iret) grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) vector(1)=.false. call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector) do kk=1,grd%iglobal i=grd%ltosi(kk) j=grd%ltosj(kk) grid3(i,j,1)=work1(kk)-grid3(i,j,1) end do call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector) do j=1,latb do i=1,lonb grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1) end do end do rwork1d = reshape(grid_b,(/size(rwork1d)/)) else call load_grid(work1,grid) rwork1d = reshape(grid,(/size(rwork1d)/)) end if if (use_fv3_aero) then call nemsio_writerecv(gfileo,'bc2','mid layer',k,rwork1d,iret=iret) else rwork1d=rwork1d*1.0e-9_r_kind ! convert ug/kg back to kg/kg call nemsio_writerecv(gfileo,'bcphilic','mid layer',k,rwork1d,iret=iret) end if if (iret /= 0) call error_msg(trim(my_name),trim(filename),'bcphilic','write',istop,iret) endif end do ! output extra variables if FV3-Chem if (use_fv3_aero) then do k=1,grd%nsig if (mype == mype_out) then call nemsio_readrecv(gfile,'pp25','mid layer',k,rwork1d,iret=iret) call nemsio_writerecv(gfileo,'pp25','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pp25','write',istop,iret) end if end do do k=1,grd%nsig if (mype == mype_out) then call nemsio_readrecv(gfile,'pp10','mid layer',k,rwork1d,iret=iret) call nemsio_writerecv(gfileo,'pp10','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pp10','write',istop,iret) end if end do end if end if ! end if laeroana_gocart ! ! Deallocate local array ! if (mype==mype_out) then if (diff_res .or. lupp .or. imp_physics == 11) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid3,clons,slons) if (imp_physics == 11) deallocate(grid3b) call nemsio_close(gfile,iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),null,'close',istop,iret) call nemsio_close(gfileo,iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret) ! ! Deallocate local array ! deallocate(rwork1d,rwork1d1) ! write(6,'(a,'': atm anal written for lonb,latb,levs= '',3i6,'',valid hour= '',f4.1,'',idate= '',4i5)') & trim(my_name),lonb,latb,levs,fhour,odate endif end subroutine write_atm_ subroutine write_sfc_ (filename,mype_sfc,dsfct) !$$$ subprogram documentation block ! . . . ! subprogram: write_nemssfc --- Write surface analysis to file ! ! prgmmr: Huang org: np23 date: 2010-02-22 ! ! abstract: This routine writes the updated surface analysis. At ! this point (20101020) the only surface field update by ! the gsi is the skin temperature. The current (20101020) ! GDAS setup does use the updated surface file. Rather, ! the output from surface cycle is used as the surface ! analysis for subsequent NEMS/GFS runs. ! ! The routine gathers surface fields from subdomains, ! reformats the data records, and then writes each record ! to the output file. ! ! Since the gsi only update the skin temperature, all ! other surface fields are simply read from the guess ! surface file and written to the analysis file. ! ! program history log: ! 2010-02-22 Huang Initial version. Based on write_gfssfc ! 2011-04-01 Huang change type of buffer2, grid2 from single to r_kind ! 2013-10-25 todling - reposition ltosi and others to commvars ! ! input argument list: ! filename - file to open and write to ! dsfct - delta skin temperature ! mype_sfc - mpi task to write output file ! ! output argument list: ! ! attributes: ! language: f90 ! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP ! !$$$ end documentation block ! !USES: use kinds, only: r_kind,i_kind,r_single use mpimod, only: mpi_rtype use mpimod, only: mpi_comm_world use mpimod, only: ierror use mpimod, only: mype use gridmod, only: nlat,nlon use gridmod, only: lat1,lon1 use gridmod, only: lat2,lon2 use gridmod, only: iglobal use gridmod, only: ijn use gridmod, only: displs_g use gridmod, only: itotsub use gridmod, only: rlats,rlons,rlats_sfc,rlons_sfc use general_commvars_mod, only: ltosi,ltosj use obsmod, only: iadate use constants, only: zero use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close,nemsio_readrecv use nemsio_module, only: nemsio_gfile,nemsio_getfilehead use nemsio_module, only: nemsio_readrec, nemsio_writerec, nemsio_writerecv implicit none ! !INPUT PARAMETERS: character(24) ,intent(in ) :: filename ! file to open and write to real(r_kind),dimension(lat2,lon2),intent(in ) :: dsfct ! delta skin temperature integer(i_kind) ,intent(in ) :: mype_sfc ! mpi task to write output file ! !OUTPUT PARAMETERS: !------------------------------------------------------------------------- ! Declare local parameters character( 6),parameter:: fname_ges='sfcf06' ! Declare local variables character(len=120) :: my_name = 'WRITE_NEMSSFC' character(len=1) :: null = ' ' integer(i_kind),dimension(7):: idate, jdate integer(i_kind),dimension(4):: odate integer(i_kind) :: i, j, ip1, jp1, ilat, ilon, jj, mm1 integer(i_kind) :: nlatm2, n, nrec, lonb, latb, iret integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd integer(i_kind) :: istop = 105 real(r_kind) :: fhour real(r_kind),dimension(lat1,lon1):: sfcsub real(r_kind),dimension(nlon,nlat):: grid real(r_kind),dimension(max(iglobal,itotsub)):: sfcall real(r_kind),allocatable,dimension(:,:) :: tsea real(r_kind),allocatable,dimension(:) :: rwork1d real(r_single),dimension(nlon,nlat):: buffer real(r_single),allocatable,dimension(:,:) :: buffer2,grid2 type(nemsio_gfile) :: gfile, gfileo !***************************************************************************** ! Initialize local variables mm1=mype+1 nlatm2=nlat-2 ! Gather skin temperature information from all tasks. do j=1,lon1 jp1 = j+1 do i=1,lat1 ip1 = i+1 sfcsub(i,j)=dsfct(ip1,jp1) end do end do call mpi_gatherv(sfcsub,ijn(mm1),mpi_rtype,& sfcall,ijn,displs_g,mpi_rtype,mype_sfc,& mpi_comm_world,ierror) ! Only MPI task mype_sfc writes the surface file. if (mype==mype_sfc) then ! Reorder updated skin temperature to output format do i=1,iglobal ilon=ltosj(i) ilat=ltosi(i) grid(ilon,ilat)=sfcall(i) end do do j=1,nlat jj=nlat-j+1 do i=1,nlon buffer(i,j)=grid(i,jj) end do end do ! Read surface guess file call nemsio_init(iret=iret) if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret) call nemsio_open(gfile,fname_ges,'read',iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),null,'open',istop,iret) ! call nemsio_getfilehead(gfile, nrec=nrec, idate=idate, dimx=lonb, & dimy=latb, nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, & nfsecondd=nfsecondd, iret=iret) ! ! Replace header record date with analysis time from iadate ! jdate(1) = iadate(1) ! analysis year jdate(2) = iadate(2) ! analysis month jdate(3) = iadate(3) ! analysis day jdate(4) = iadate(4) ! analysis hour jdate(5) = iadate(5) ! analysis minute jdate(5) = 0 ! analysis minute jdate(6) = 0 ! analysis scaled seconds jdate(7) = idate(7) ! analysis seconds multiplier nfhour=0 ! new forecast hour, zero at analysis time nfminute=0 nfsecondn=0 nfsecondd=100 ! default for denominator fhour = zero odate(1) = jdate(4) !hour odate(2) = jdate(2) !month odate(3) = jdate(3) !day odate(4) = jdate(1) !year ! ! Start to write output sfc file : filename ! open new output file with new header gfileo with "write" access. ! Use this call to update header as well ! ! gfileo=gfile ! copy input header info to output header info ! need to do this before nemsio_close(gfile) call nemsio_open(gfileo,filename,'write',iret=iret, idate=jdate, nfhour=nfhour,& nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd ) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret) ! ! First copy entire data from fname_ges to filename, then do selective update ! allocate(rwork1d(lonb*latb)) allocate(buffer2(lonb,latb)) allocate(grid2(lonb,latb)) allocate(tsea(lonb,latb)) do n = 1, nrec call nemsio_readrec (gfile, n,rwork1d,iret=iret) if ( iret /= 0 ) write(6,*) 'readrec nrec = ', n, ' Status = ', iret call nemsio_writerec(gfileo,n,rwork1d,iret=iret) if ( iret /= 0 ) write(6,*) 'writerec nrec = ', n, ' Status = ', iret end do ! ! Only sea surface temperature will be updated in the SFC files ! call nemsio_readrecv(gfile,'tmp','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),'tmp','read',istop,iret) tsea=reshape(rwork1d,(/size(tsea,1),size(tsea,2)/)) if ( (latb /= nlatm2) .or. (lonb /= nlon) ) then write(6,*)trim(my_name),': different grid dimensions analysis', & ' vs sfc. interpolating sfc temperature nlon,nlat-2=',nlon, & nlatm2,' -vs- sfc file lonb,latb=',lonb,latb call intrp22(buffer, rlons,rlats,nlon,nlat, & buffer2,rlons_sfc,rlats_sfc,lonb,latb) else do j=1,latb do i=1,lonb buffer2(i,j)=buffer(i,j+1) end do end do endif grid2 = tsea + buffer2 rwork1d = reshape( grid2,(/size(rwork1d)/) ) deallocate(buffer2) ! update tsea record call nemsio_writerecv(gfileo,'tmp','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','write',istop,iret) deallocate(rwork1d) call nemsio_close(gfile, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),null,'close',istop,iret) call nemsio_close(gfileo,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret) write(6,'(a,'': sfc anal written for lonb,latb= '',2i6,'',valid hour= '',f4.1,'',idate= '',4i5)') & trim(my_name),lonb,latb,fhour,odate endif end subroutine write_sfc_ subroutine write_sfc_nst_ (mype_so,dsfct) !$$$ subprogram documentation block ! . . . ! subprogram: write_sfc_nst --- Write both sfc and nst surface analysis to file ! ! prgmmr: Huang org: np23 date: 2011-11-01 ! ! abstract: This routine writes the sfc & nst analysis files and is nst_gsi dependent. ! Tr (foundation temperature), instead of skin temperature, is the analysis variable. ! nst_gsi > 2: Tr analysis is on ! nst_gsi <= 2: Tr analysis is off ! ! The routine gathers Tr field from subdomains, ! reformats the data records, and then writes each record ! to the output files. ! ! Since the gsi only update the Tr temperature, all ! other fields in surface are simply read from the guess ! files and written to the analysis file. ! ! program history log: ! 2011-11-01 Huang initial version based on routine write_gfs_sfc_nst ! 2013-10-25 todling - reposition ltosi and others to commvars ! 2016-01-01 li - update write_sfc_nst_ (nemsio) as for write_gfs_sfc_nst (sfcio) ! ! input argument list: ! dsfct - delta skin temperature ! mype_so - mpi task to write output file ! ! output argument list: ! ! attributes: ! language: f90 ! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP ! !$$$ end documentation block ! !USES: use kinds, only: r_kind,i_kind,r_single use mpimod, only: mpi_rtype,mpi_itype use mpimod, only: mpi_comm_world use mpimod, only: ierror use mpimod, only: mype use gridmod, only: nlat,nlon use gridmod, only: lat1,lon1 use gridmod, only: lat2,lon2 use gridmod, only: nlat_sfc,nlon_sfc use gridmod, only: iglobal use gridmod, only: ijn use gridmod, only: displs_g use gridmod, only: itotsub use general_commvars_mod, only: ltosi,ltosj use obsmod, only: iadate use constants, only: zero,two,tfrozen,z_w_max use constants, only: zero_single use guess_grids, only: isli2 use gsi_nstcouplermod, only: nst_gsi,zsea1,zsea2 use gridmod, only: rlats,rlons,rlats_sfc,rlons_sfc use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close,nemsio_readrecv use nemsio_module, only: nemsio_gfile,nemsio_getfilehead use nemsio_module, only: nemsio_readrec, nemsio_writerec, nemsio_writerecv implicit none ! !INPUT PARAMETERS: real(r_kind),dimension(lat2,lon2),intent(in ) :: dsfct ! delta skin temperature integer(i_kind) ,intent(in ) :: mype_so ! mpi task to write output file ! !OUTPUT PARAMETERS: !------------------------------------------------------------------------- ! Declare local parameters character(6), parameter:: fname_sfcges = 'sfcf06' character(6), parameter:: fname_sfcgcy = 'sfcgcy' character(6), parameter:: fname_sfctsk = 'sfctsk' character(6), parameter:: fname_sfcanl = 'sfcanl' character(6), parameter:: fname_nstges = 'nstf06' character(6), parameter:: fname_nstanl = 'nstanl' character(6), parameter:: fname_dtfanl = 'dtfanl' ! Declare local variables integer(i_kind), parameter:: io_dtfanl = 54 integer(i_kind), parameter:: nprep=15 real(r_kind),parameter :: houra = zero_single character(len=120) :: my_name = 'WRITE_SFC_NST' character(len=1) :: null = ' ' integer(i_kind),dimension(7):: idate, jdate integer(i_kind),dimension(4):: odate integer(i_kind) :: i, j, ip1, jp1, ilat, ilon, mm1 integer(i_kind) :: lonb, latb, nlatm2, n, nrec_sfc, nrec_nst, iret integer(i_kind) :: lonb_nst, latb_nst integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd integer(i_kind) :: istop = 106 real(r_kind) :: fhour real(r_single) :: r_zsea1,r_zsea2 real(r_kind), dimension(lat1,lon1):: dsfct_sub integer(i_kind), dimension(lat1,lon1):: isli_sub real(r_kind), dimension(max(iglobal,itotsub)):: dsfct_all integer(i_kind), dimension(max(iglobal,itotsub)):: isli_all real(r_kind), dimension(nlat,nlon):: dsfct_glb,dsfct_tmp integer(i_kind), dimension(nlat,nlon):: isli_glb,isli_tmp real(r_kind), dimension(nlat_sfc,nlon_sfc) :: dsfct_gsi integer(i_kind), dimension(nlat_sfc,nlon_sfc) :: isli_gsi real(r_kind), dimension(nlon_sfc,nlat_sfc-2):: dsfct_anl real(r_single), dimension(nlon_sfc,nlat_sfc-2):: dtzm real(r_single), dimension(nlat_sfc,nlon_sfc) :: work real(r_single), allocatable, dimension(:,:) :: tsea,xt,xs,xu,xv,xz,zm,xtts,xzts,dt_cool,z_c, & c_0,c_d,w_0,w_d,d_conv,ifd,tref,qrain real(r_single), allocatable, dimension(:,:) :: slmsk_ges,slmsk_anl real(r_single), allocatable, dimension(:) :: rwork1d type(nemsio_gfile) :: gfile_sfcges,gfile_sfcgcy,gfile_nstges,gfile_sfctsk,gfile_sfcanl,gfile_nstanl !***************************************************************************** ! Initialize local variables mm1=mype+1 nlatm2=nlat-2 ! ! Extract the analysis increment and surface mask in subdomain without the buffer ! do j=1,lon1 jp1 = j+1 do i=1,lat1 ip1 = i+1 dsfct_sub(i,j) = dsfct(ip1,jp1) isli_sub (i,j) = isli2(ip1,jp1) end do end do ! ! Gather global analysis increment and surface mask info from subdomains ! call mpi_gatherv(dsfct_sub,ijn(mm1),mpi_rtype,& dsfct_all,ijn,displs_g,mpi_rtype,mype_so ,& mpi_comm_world,ierror) call mpi_gatherv(isli_sub,ijn(mm1),mpi_itype,& isli_all,ijn,displs_g,mpi_itype,mype_so ,& mpi_comm_world,ierror) ! Only MPI task mype_so writes the surface file. if (mype==mype_so ) then write(*,'(a,5(1x,a6))') 'write_nems_sfc_nst:',fname_sfcges,fname_nstges,fname_sfctsk,fname_sfcanl,fname_nstanl ! ! get Tf analysis increment and surface mask at analysis (lower resolution) grids ! do i=1,iglobal ilon=ltosj(i) ilat=ltosi(i) dsfct_glb(ilat,ilon) = dsfct_all(i) isli_glb (ilat,ilon) = isli_all (i) end do ! ! write dsfct_anl to a data file for later use (at eupd step at present) ! open(io_dtfanl,file=fname_dtfanl,form='unformatted') write(io_dtfanl) nlon,nlat write(io_dtfanl) dsfct_glb write(io_dtfanl) isli_glb ! Initiate nemsio call nemsio_init(iret=iret) if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret) ! open nsst guess file call nemsio_open(gfile_nstges,trim(fname_nstges),'read',iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),null,'open',istop,iret) ! open surface guess file call nemsio_open(gfile_sfcges,trim(fname_sfcges),'read',iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcges),null,'open',istop,iret) ! open surface gcycle file call nemsio_open(gfile_sfcgcy,trim(fname_sfcgcy),'read',iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcgcy),null,'open',istop,iret) ! read a few surface guess file header records call nemsio_getfilehead(gfile_sfcges, nrec=nrec_sfc, idate=idate, & dimx=lonb, dimy=latb, nfhour=nfhour, nfminute=nfminute, & nfsecondn=nfsecondn, nfsecondd=nfsecondd, iret=iret) ! read some nsst guess file header records (dimensions) call nemsio_getfilehead(gfile_nstges, nrec=nrec_nst, dimx=lonb_nst,dimy=latb_nst,iret=iret) write(6,*) 'nrec_sfc, nrec_nst = ',nrec_sfc, nrec_nst ! check the dimensions consistency in sfc, nst files and the used. if ( latb /= latb_nst .or. lonb /= lonb_nst ) then write(6,*) 'Inconsistent dimension for sfc & nst files. latb,lonb : ',latb,lonb, & 'latb_nst,lonb_nst : ',latb_nst,lonb_nst call stop2(80) endif if ( nlat_sfc /= latb+2 .or. nlon_sfc /= lonb ) then write(6,*) 'Inconsistent dimension for used and read. nlat_sfc,nlon_sfc : ',nlat_sfc,nlon_sfc, & 'latb+2,lonb :',latb+2,lonb call stop2(81) endif ! allocate(slmsk_ges(lonb,latb),slmsk_anl(lonb,latb)) allocate(rwork1d(lonb*latb)) ! read slmsk in fname_sfcges to get slmsk_ges call nemsio_readrecv(gfile_sfcges, 'land', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcges),'land','read',istop,iret) slmsk_ges=reshape(rwork1d,(/size(slmsk_ges,1),size(slmsk_ges,2)/)) ! read slmsk in fname_sfcgcy to get slmsk_anl call nemsio_readrecv(gfile_sfcgcy, 'land', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcgcy),'land','read',istop,iret) slmsk_anl=reshape(rwork1d,(/size(slmsk_anl,1),size(slmsk_anl,2)/)) ! ! Replace header record date with analysis time from iadate ! jdate(1) = iadate(1) ! analysis year jdate(2) = iadate(2) ! analysis month jdate(3) = iadate(3) ! analysis day jdate(4) = iadate(4) ! analysis hour jdate(5) = iadate(5) ! analysis minute jdate(5) = 0 ! analysis minute jdate(6) = 0 ! analysis scaled seconds jdate(7) = idate(7) ! analysis seconds multiplier nfhour=0 ! new forecast hour, zero at analysis time nfminute=0 nfsecondn=0 nfsecondd=100 ! default for denominator fhour = zero odate(1) = jdate(4) !hour odate(2) = jdate(2) !month odate(3) = jdate(3) !day odate(4) = jdate(1) !year if ( (latb /= nlatm2) .or. (lonb /= nlon) ) then write(6,*)'WRITE_NEMSIO_SFC_NST: different grid dimensions analysis vs sfc. interpolating sfc temperature ',& ', nlon,nlat-2=',nlon,nlatm2,' -vs- sfc file lonb,latb=',lonb,latb write(6,*) ' WRITE_NEMSIO_SFC_NST, nlon_sfc,nlat_sfc : ', nlon_sfc,nlat_sfc ! ! Get the expanded values for a surface type (0 = water now) and the new mask ! call int2_msk_glb_prep(dsfct_glb,isli_glb,dsfct_tmp,isli_tmp,nlat,nlon,0,nprep) ! ! Get updated/analysis surface mask info from sfcgcy file ! call tran_gfssfc(slmsk_anl,work,lonb,latb) do j=1,lonb do i=1,latb+2 isli_gsi(i,j) = nint(work(i,j)) end do end do ! ! Interpolate dsfct_tmp(nlat,nlon) to dsfct_gsi(nlat_sfc,nlon_sfc) with surface mask accounted ! call int22_msk_glb(dsfct_tmp,isli_tmp,rlats,rlons,nlat,nlon, & dsfct_gsi,isli_gsi,rlats_sfc,rlons_sfc,nlat_sfc,nlon_sfc,0) ! ! transform the dsfct_gsi(latb+2,lonb) to dsfct_anl(lonb,latb) for sfc file format ! do j = 1, latb do i = 1, lonb dsfct_anl(i,j) = dsfct_gsi(latb+2-j,i) end do end do else ! ! transform the dsfct_glb(nlat,nlon) to dsfct_anl(lonb,latb) for sfc file ! format when nlat == latb-2 & nlon = lonb ! do j=1,latb do i=1,lonb dsfct_anl(i,j)=dsfct_glb(latb+1-j,i) end do end do endif ! if ( (latb /= nlatm2) .or. (lonb /= nlon) ) then ! ! Start to write output sfc file : fname_sfcanl & fname_nstanl ! open new output file with new header gfile_sfcanl and gfile_nstanl with "write" access. ! Use this call to update header as well ! ! copy input header info to output header info for sfcanl, need to do this before nemsio_close(gfile) ! gfile_sfcanl=gfile_sfcgcy ! open nemsio sfcanl call nemsio_open(gfile_sfcanl,trim(fname_sfcanl),'write',iret=iret, idate=jdate, nfhour=nfhour,& nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd ) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcanl),null,'open',istop,iret) gfile_sfctsk=gfile_sfcgcy ! open nemsio sfctsk call nemsio_open(gfile_sfctsk,trim(fname_sfctsk),'write',iret=iret, idate=jdate, nfhour=nfhour,& nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd ) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfctsk),null,'open',istop,iret) ! ! copy input header info to output header info for nstanl, need to do this before nemsio_close(gfile) ! gfile_nstanl=gfile_nstges ! open nemsio nstanl call nemsio_open(gfile_nstanl,trim(fname_nstanl),'write',iret=iret, idate=jdate, nfhour=nfhour,& nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd ) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),null,'open',istop,iret) ! Allocate work array (rwork1d) and tsea in sfc file allocate(tsea(lonb,latb)) ! Allocate nsst variables allocate(xt(lonb,latb)) allocate(xs(lonb,latb)) allocate(xu(lonb,latb)) allocate(xv(lonb,latb)) allocate(xz(lonb,latb)) allocate(zm(lonb,latb)) allocate(xtts(lonb,latb)) allocate(xzts(lonb,latb)) allocate(dt_cool(lonb,latb)) allocate(z_c(lonb,latb)) allocate(c_0(lonb,latb)) allocate(c_d(lonb,latb)) allocate(w_0(lonb,latb)) allocate(w_d(lonb,latb)) allocate(d_conv(lonb,latb)) allocate(ifd(lonb,latb)) allocate(tref(lonb,latb)) allocate(qrain(lonb,latb)) ! ! First copy entire data from sfcgcy to fname_anl, then do selective update ! ! read the nrec_sfc variables from sfcgcy and then write then to sfcanl ! do n = 1, nrec_sfc call nemsio_readrec(gfile_sfcgcy,n,rwork1d,iret=iret) if ( iret /= 0 ) write(6,*) 'readrec for gfile_sfcgcy, nrec_sfc = ', n, ' Status = ', iret call nemsio_writerec(gfile_sfcanl,n,rwork1d,iret=iret) if ( iret /= 0 ) write(6,*) 'writerec for gfile_sfcanl, nrec_sfc = ', n, ' Status = ', iret call nemsio_writerec(gfile_sfctsk,n,rwork1d,iret=iret) if ( iret /= 0 ) write(6,*) 'writerec for gfile_sfctsk, nrec_sfc = ', n, ' Status = ', iret end do write(*,*) 'read gfile_sfcgcy, and the write to gfile_sfcanl, gfile_sfctsk' ! ! For sfcanl, Only tsea (sea surface temperature) will be updated in the SFC ! Need values from nstges for tref update ! read tsea from sfcges call nemsio_readrecv(gfile_sfcges,'tmp','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcges),'tmp','read',istop,iret) tsea=reshape(rwork1d,(/size(tsea,1),size(tsea,2)/)) ! For nstanl, Only tref (foundation temperature) is updated by analysis ! others are updated for snow melting case ! read 18 nsst variables from nstges ! xt call nemsio_readrecv(gfile_nstges, 'xt', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xt','read',istop,iret) xt=reshape(rwork1d,(/size(xt,1),size(xt,2)/)) ! xs call nemsio_readrecv(gfile_nstges, 'xs', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xs','read',istop,iret) xs=reshape(rwork1d,(/size(xs,1),size(xs,2)/)) ! xu call nemsio_readrecv(gfile_nstges, 'xu', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xu','read',istop,iret) xu=reshape(rwork1d,(/size(xu,1),size(xu,2)/)) ! xv call nemsio_readrecv(gfile_nstges, 'xv', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xv','read',istop,iret) xv=reshape(rwork1d,(/size(xv,1),size(xv,2)/)) ! xz call nemsio_readrecv(gfile_nstges, 'xz', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xz','read',istop,iret) xz=reshape(rwork1d,(/size(xz,1),size(xz,2)/)) ! zm call nemsio_readrecv(gfile_nstges, 'zm', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'zm','read',istop,iret) zm=reshape(rwork1d,(/size(zm,1),size(zm,2)/)) ! xtts call nemsio_readrecv(gfile_nstges, 'xtts', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xtts','read',istop,iret) xtts=reshape(rwork1d,(/size(xtts,1),size(xtts,2)/)) ! xzts call nemsio_readrecv(gfile_nstges, 'xzts', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xzts','read',istop,iret) xzts=reshape(rwork1d,(/size(xzts,1),size(xzts,2)/)) ! dt_cool call nemsio_readrecv(gfile_nstges, 'dtcool','sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'dt_cool','read',istop,iret) dt_cool=reshape(rwork1d,(/size(dt_cool,1),size(dt_cool,2)/)) ! z_c call nemsio_readrecv(gfile_nstges, 'zc','sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'zc','read',istop,iret) z_c=reshape(rwork1d,(/size(z_c,1),size(z_c,2)/)) ! c_0 call nemsio_readrecv(gfile_nstges, 'c0','sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'c0','read',istop,iret) c_0=reshape(rwork1d,(/size(c_0,1),size(c_0,2)/)) ! c_d call nemsio_readrecv(gfile_nstges, 'cd','sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'cd','read',istop,iret) c_d=reshape(rwork1d,(/size(c_d,1),size(c_d,2)/)) ! w_0 call nemsio_readrecv(gfile_nstges, 'w0','sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'w0','read',istop,iret) w_0=reshape(rwork1d,(/size(w_0,1),size(w_0,2)/)) ! w_d call nemsio_readrecv(gfile_nstges, 'wd','sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'wd','read',istop,iret) w_d=reshape(rwork1d,(/size(w_d,1),size(w_d,2)/)) ! tref call nemsio_readrecv(gfile_nstges, 'tref', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'tref','read',istop,iret) tref=reshape(rwork1d,(/size(tref,1),size(tref,2)/)) ! d_conv call nemsio_readrecv(gfile_nstges, 'dconv', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'dconv','read',istop,iret) d_conv=reshape(rwork1d,(/size(d_conv,1),size(d_conv,2)/)) ! ifd call nemsio_readrecv(gfile_nstges, 'ifd', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'ifd','read',istop,iret) ifd=reshape(rwork1d,(/size(ifd,1),size(ifd,2)/)) ! qrain call nemsio_readrecv(gfile_nstges, 'qrain', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'qrain','read',istop,iret) qrain=reshape(rwork1d,(/size(qrain,1),size(qrain,2)/)) ! ! update tref (in nst file) & tsea (in the surface file) when Tr analysis is on ! reset NSSTM variables for new open water grids ! if ( nst_gsi > 2 ) then ! ! For the new open water (sea ice just melted) grids, (1) set dsfct_anl = zero; (2) reset the NSSTM variables ! ! Notes: slmsk_ges is the mask of the background ! slmsk_anl is the mask of the analysis ! where ( slmsk_anl(:,:) == zero .and. slmsk_ges(:,:) == two ) dsfct_anl(:,:) = zero xt(:,:) = zero xs(:,:) = zero xu(:,:) = zero xv(:,:) = zero xz(:,:) = z_w_max zm(:,:) = zero xtts(:,:) = zero xzts(:,:) = zero dt_cool(:,:) = zero z_c(:,:) = zero c_0(:,:) = zero c_d(:,:) = zero w_0(:,:) = zero w_d(:,:) = zero d_conv(:,:) = zero ifd(:,:) = zero tref(:,:) = tfrozen qrain(:,:) = zero end where ! ! update analysis variable: Tref (foundation temperature) for nst file ! where ( slmsk_anl(:,:) == zero ) tref(:,:) = max(tref(:,:) + dsfct_anl(:,:),tfrozen) elsewhere tref(:,:) = tsea(:,:) end where ! ! update SST: tsea for sfc file with NSST profile ! r_zsea1 = 0.001_r_single*real(zsea1) r_zsea2 = 0.001_r_single*real(zsea2) call dtzm_2d(xt,xz,dt_cool,z_c,slmsk_anl,r_zsea1,r_zsea2,lonb,latb,dtzm) where ( slmsk_anl(:,:) == zero ) tsea(:,:) = max(tref(:,:) + dtzm(:,:), tfrozen) end where else ! when (nst_gsi <= 2) do j=1,latb do i=1,lonb tref(i,j) = tsea(i,j) ! keep tref as tsea before analysis end do end do ! ! For the new open water (sea ice just melted) grids, reset the NSSTM variables ! where ( slmsk_anl(:,:) == zero .and. slmsk_ges(:,:) == two ) xt(:,:) = zero xs(:,:) = zero xu(:,:) = zero xv(:,:) = zero xz(:,:) = z_w_max zm(:,:) = zero xtts(:,:) = zero xzts(:,:) = zero dt_cool(:,:) = zero z_c(:,:) = zero c_0(:,:) = zero c_d(:,:) = zero w_0(:,:) = zero w_d(:,:) = zero d_conv(:,:) = zero ifd(:,:) = zero tref(:,:) = tfrozen qrain(:,:) = zero end where ! ! update tsea when NO Tf analysis ! do j=1,latb do i=1,lonb tsea(i,j) = max(tsea(i,j) + dsfct_anl(i,j),tfrozen) end do end do endif ! if ( nst_gsi > 2 ) then ! ! update tsea record in sfcanl ! rwork1d = reshape(tsea, (/size(rwork1d)/) ) call nemsio_writerecv(gfile_sfcanl,'tmp','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcanl),'tmp','write',istop,iret) write(6,100) fname_sfcanl,lonb,latb,houra,iadate(1:4),iret 100 format(' WRITE_NEMSIO_SFC_NST: update tsea in ',a6,2i6,1x,f4.1,4(i4,1x),' with iret=',i2) ! ! update tsea record in sfctsk ! rwork1d = reshape(tsea, (/size(rwork1d)/) ) call nemsio_writerecv(gfile_sfctsk,'tmp','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfctsk),'tmp','write',istop,iret) write(6,101) fname_sfctsk,lonb,latb,houra,iadate(1:4),iret 101 format(' WRITE_NEMSIO_SFC_NST: update tsea in ',a6,2i6,1x,f4.1,4(i4,1x),' with iret=',i2) ! ! update nsst records in nstanl ! ! slmsk rwork1d = reshape( slmsk_anl,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'land','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'land','write',istop,iret) ! xt rwork1d = reshape( xt,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'xt','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xt','write',istop,iret) ! xs rwork1d = reshape( xs,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'xs','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xs','write',istop,iret) ! xu rwork1d = reshape( xu,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'xu','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xu','write',istop,iret) ! xv rwork1d = reshape( xv,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'xv','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xv','write',istop,iret) ! xz rwork1d = reshape( xz,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'xz','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xz','write',istop,iret) ! zm rwork1d = reshape( zm,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'zm','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'zm','write',istop,iret) ! xtts rwork1d = reshape( xtts,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'xtts','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xtts','write',istop,iret) ! xzts rwork1d = reshape( xzts,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'xzts','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xzts','write',istop,iret) ! z_0 rwork1d = reshape( dt_cool,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'dtcool','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'dtcool','write',istop,iret) ! z_c rwork1d = reshape( z_c,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'zc','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'zc','write',istop,iret) ! c_0 rwork1d = reshape( c_0,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'c0','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'c0','write',istop,iret) ! c_d rwork1d = reshape( c_d,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'cd','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'cd','write',istop,iret) ! w_0 rwork1d = reshape( w_0,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'w0','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'w0','write',istop,iret) ! w_d rwork1d = reshape( w_d,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'wd','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'wd','write',istop,iret) ! d_conv rwork1d = reshape( d_conv,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'dconv','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'dconv','write',istop,iret) ! ifd rwork1d = reshape( ifd,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'ifd','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'ifd','write',istop,iret) ! tref rwork1d = reshape( tref,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'tref','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'tref','write',istop,iret) ! qrain rwork1d = reshape( qrain,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'qrain','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'qrain','write',istop,iret) write(6,200) fname_nstanl,lonb,latb,houra,iadate(1:4),iret 200 format(' WRITE_NEMSIO_SFC_NST: update variables in ',a6,2i6,1x,f4.1,4(i4,1x),' with iret=',i2) deallocate(xt,xs,xu,xv,xz,zm,xtts,xzts,dt_cool,z_c,c_0,c_d,w_0,w_d,d_conv,ifd,tref,qrain) deallocate(rwork1d) call nemsio_close(gfile_sfcges, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcges),null,'close',istop,iret) call nemsio_close(gfile_sfcgcy, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcgcy),null,'close',istop,iret) call nemsio_close(gfile_nstges, iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),null,'close',istop,iret) call nemsio_close(gfile_sfcanl,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcanl),null,'close',istop,iret) call nemsio_close(gfile_nstanl,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),null,'close',istop,iret) call nemsio_close(gfile_sfctsk,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfctsk),null,'close',istop,iret) write(6,'(a,'': nemsio sfc_nst anal written for lonb,latb= '',2i6,'',valid hour= '',f4.1,'',idate= '',4i5)') & trim(my_name),lonb,latb,fhour,odate endif end subroutine write_sfc_nst_ subroutine error_msg_(sub_name,file_name,var_name,action,stop_code,error_code,lprint) use mpimod, only: mype use kinds, only: i_kind implicit none character(len=*), intent(in) :: sub_name,file_name,var_name,action integer(i_kind), intent(in) :: stop_code, error_code logical, optional,intent(in) :: lprint if ( mype == 0 .or. present(lprint) ) then select case (trim(action)) case('init') write(6,'(a,'': PROBLEM with nemsio_init, Status = '', i3)') & trim(sub_name), error_code case('open') write(6,'(a,'': ***ERROR*** problem opening file '',a,'', Status = '', i3)') & trim(sub_name), trim(file_name), error_code case('close') write(6,'(a,'': ***ERROR*** problem closing file '',a,'', Status = '', i3)') & trim(sub_name), trim(file_name), error_code case default write(6,'(a,'': ***ERROR*** '',a,tr1,a,'',variable = '',a,'',Status = '',i3)') & trim(sub_name),trim(action),trim(file_name),trim(var_name),error_code end select end if if ( stop_code /= 0 ) call stop2(stop_code) end subroutine error_msg_ subroutine intrp22(a,rlons_a,rlats_a,nlon_a,nlat_a, & b,rlons_b,rlats_b,nlon_b,nlat_b) !$$$ subprogram documentation block ! . . . ! subprogram: intrp22 --- interpolates from one 2-d grid to another 2-d grid ! like analysis to surface grid or vice versa ! prgrmmr: li - initial version; org: np2 ! ! abstract: This routine interpolates a grid to b grid ! ! program history log: ! ! input argument list: ! rlons_a - longitudes of input array ! rlats_a - latitudes of input array ! nlon_a - number of longitude of input array ! nlat_a - number of latitude of input array ! rlons_b - longitudes of output array ! rlats_b - latitudes of output array ! nlon_b - number of longitude of output array ! nlat_b - number of latitude of output array ! a - input values ! ! output argument list: ! b - output values ! ! attributes: ! language: f90 ! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP ! !$$$ end documentation block ! !USES: use kinds, only: r_kind,i_kind,r_single use constants, only: zero,one implicit none ! !INPUT PARAMETERS: integer(i_kind) ,intent(in ) :: nlon_a,nlat_a,nlon_b,nlat_b real(r_kind), dimension(nlon_a) ,intent(in ) :: rlons_a real(r_kind), dimension(nlat_a) ,intent(in ) :: rlats_a real(r_kind), dimension(nlon_b) ,intent(in ) :: rlons_b real(r_kind), dimension(nlat_b) ,intent(in ) :: rlats_b real(r_single), dimension(nlon_a,nlat_a),intent(in ) :: a ! !OUTPUT PARAMETERS: real(r_single), dimension(nlon_b,nlat_b),intent( out) :: b ! Declare local variables integer(i_kind) i,j,ix,iy,ixp,iyp real(r_kind) dx1,dy1,dx,dy,w00,w01,w10,w11,bout,dlat,dlon !***************************************************************************** b=zero ! Loop over all points to get interpolated value do j=1,nlat_b dlat=rlats_b(j) call grdcrd1(dlat,rlats_a,nlat_a,1) iy=int(dlat) iy=min(max(1,iy),nlat_a) dy =dlat-iy dy1 =one-dy iyp=min(nlat_a,iy+1) do i=1,nlon_b dlon=rlons_b(i) call grdcrd1(dlon,rlons_a,nlon_a,1) ix=int(dlon) dx =dlon-ix dx=max(zero,min(dx,one)) dx1 =one-dx w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy ix=min(max(0,ix),nlon_a) ixp=ix+1 if(ix==0) ix=nlon_a if(ixp==nlon_a+1) ixp=1 bout=w00*a(ix,iy)+w01*a(ix,iyp)+w10*a(ixp,iy)+w11*a(ixp,iyp) b(i,j)=bout end do end do ! End of routine return end subroutine intrp22 subroutine tran_gfssfc(ain,aout,lonb,latb) !$$$ subprogram documentation block ! . . . . ! subprogram: tran_gfssfc transform gfs surface file to analysis grid ! prgmmr: derber org: np2 date: 2003-04-10 ! ! abstract: transform gfs surface file to analysis grid ! ! program history log: ! 2012-31-38 derber - initial routine ! ! input argument list: ! ain - input surface record on processor iope ! lonb - input number of longitudes ! latb - input number of latitudes ! ! output argument list: ! aout - output transposed surface record ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind,r_single use constants, only: zero use sfcio_module, only: sfcio_realkind implicit none ! Declare passed variables integer(i_kind) ,intent(in ) :: lonb,latb real(sfcio_realkind),dimension(lonb,latb),intent(in ) :: ain real(r_single),dimension(latb+2,lonb),intent(out) :: aout ! Declare local variables integer(i_kind) i,j real(r_kind) sumn,sums ! of surface guess array sumn = zero sums = zero do i=1,lonb sumn = ain(i,1) + sumn sums = ain(i,latb) + sums end do sumn = sumn/float(lonb) sums = sums/float(lonb) ! Transfer from local work array to surface guess array do j = 1,lonb aout(1,j)=sums do i=2,latb+1 aout(i,j) = ain(j,latb+2-i) end do aout(latb+2,j)=sumn end do return end subroutine tran_gfssfc end module ncepnems_io