program getsfcnstensupdp !$$$ main program documentation block ! ! program: getsfcnstenupdp update sfc and nst files for ensemble ! ! prgmmr: Xu Li org: EMC date: 2014-05-01 ! ! abstract: update sfc & nst file for ensemble ! ! program history log: ! 2014-05-01 Initial version. ! 2016-02-15 Add to read nemsio ! 2016-08-18 Fix two bugs (tsea, dtf_ens) ! 2016-11-18 change nst mask name from slmsk to land ! ! usage: ! input files: ! ! output files: ! ! attributes: ! language: f95 ! !$$$ use mpi use kinds, only: r_kind,i_kind,r_single use constants, only: two,half,zero,z_w_max,tfrozen,init_constants_derived,pi use sfcio_module, only: sfcio_srohdc,sfcio_head,sfcio_data,sfcio_swohdc use nstio_module, only: nstio_srohdc,nstio_head,nstio_data,nstio_swohdc use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close use nemsio_module, only: nemsio_gfile,nemsio_getfilehead,nemsio_readrec,& nemsio_writerec,nemsio_readrecv,nemsio_writerecv implicit none logical:: nemsio, sfcio real (r_kind), parameter :: houra=zero integer(i_kind), parameter :: nprep=15 integer(i_kind) :: istop = 106 integer(i_kind), parameter :: lun_dtfanl=11,lun_nstges=21,lun_sfcges=22, & lun_sfcgcy=23,lun_nstanl=61,lun_sfcanl=62 integer(i_kind), parameter :: idrt=4 character(len=80) :: fname_dtfanl,fname_nstges,fname_sfcgcy,fname_nstanl,fname_sfcanl character(len=3) :: charnanal character(len=8) :: charbuf character(len=60) :: my_name = 'getsfcnstensupdp' character(len=1) :: null = ' ' integer(i_kind) :: mype,mype1,npe,nproc,iret integer nrec_sfc, lonb, latb, n integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd integer,dimension(7):: idate integer(i_kind) :: n_new_water,n_new_seaice integer(i_kind) :: i,j,jmax integer(i_kind) :: nanals,nst_gsi,zsea1,zsea2 integer(i_kind) :: nlon_anl,nlat_anl ! the number of lon/lat of GSI analysis grids, including two extra polar lats integer(i_kind) :: nlon_ens,nlat_ens ! the number of lon/lat of ensemble grids, including two extra polar lats integer(i_kind), allocatable, dimension(:,:) :: isli_anl,isli_epd,isli_gsi real(r_kind), allocatable, dimension(:) :: wlatx,slatx,rlats_anl,rlons_anl,rlats_ens,rlons_ens real(r_kind), allocatable, dimension(:,:) :: dtf_anl,dtf_epd,dtf_gsi,dtf_ens real(r_single), allocatable, dimension(:,:) :: work real(r_single), allocatable, dimension(:,:) :: dtzm real(r_single), allocatable, dimension(:,:) :: slmsk_ges,slmsk_ens real(r_single), allocatable, dimension(:) :: rwork1d 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) :: r_zsea1,r_zsea2 real(r_kind) :: dlon real(r_kind) :: sumn,sums type(nstio_head):: head_nst type(nstio_data):: data_nst type(sfcio_head):: head_sfcanl,head_sfcges,head_sfcgcy type(sfcio_data):: data_sfcanl,data_sfcges,data_sfcgcy type(nemsio_gfile) :: gfile_sfcges,gfile_sfcgcy,gfile_sfcanl,gfile_nstges,gfile_nstanl ! Initialize mpi ! mype is process number, npe is total number of processes. call mpi_init(iret) call MPI_Comm_rank(MPI_COMM_WORLD,mype,iret) call MPI_Comm_size(MPI_COMM_WORLD,npe,iret) call init_constants_derived if ( mype == 0 ) call w3tagb('GETSFCNSTENSUPD',2014,0921,0055,'NP25') call getarg(1,charbuf) read(charbuf,'(i3)') nanals call getarg(2,charbuf) read(charbuf,'(i1)') nst_gsi call getarg(3,charbuf) read(charbuf,'(i8)') zsea1 r_zsea1 = 0.001_r_single * real(zsea1,r_single) call getarg(4,charbuf) read(charbuf,'(i8)') zsea2 r_zsea2 = 0.001_r_single * real(zsea2,r_single) if (mype==0) then write(6,'(a)')' ' write(6,'(a)')'Command line input' write(6,'(a,i5)')' nanals = ',nanals write(6,'(a,i5)')' nst_gsi = ',nst_gsi write(6,'(a,i5)')' zsea1 = ',zsea1 write(6,'(a,i5)')' zsea2 = ',zsea2 endif ! ! read Tf analysis increment at GSI analysis grids and its grid info and surface mask ! fname_dtfanl = 'dtfanl' open(lun_dtfanl,file=trim(fname_dtfanl),form='unformatted') read(lun_dtfanl) nlon_anl,nlat_anl allocate(dtf_anl(nlat_anl,nlon_anl),isli_anl(nlat_anl,nlon_anl)) allocate(dtf_epd(nlat_anl,nlon_anl),isli_epd(nlat_anl,nlon_anl)) read(lun_dtfanl) dtf_anl read(lun_dtfanl) isli_anl ! ! determine sfcio or nemsio with sfcges_mem001 only ! sfcio=.false. nemsio=.false. call sfcio_srohdc(lun_sfcges,'sfcf06_mem001',head_sfcges,data_sfcges,iret) if ( iret == 0 ) then sfcio = .true. if (mype==0) write(6,'(3a)')'Read ','sfcf06_mem001',' in sfcio format ' lonb=head_sfcges%lonb latb=head_sfcges%latb else call nemsio_init(iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),null,null,'init',istop,iret) ! open gfile_sfcges (sfcf06_mem001) call nemsio_open(gfile_sfcges,'sfcf06_mem001','read',iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),'sfcf06_mem001',null,'open',istop,iret) if (iret == 0 ) then nemsio = .true. ! read a few surface gcycle file header records: dimensions and time parameters call nemsio_getfilehead(gfile_sfcges, nrec=nrec_sfc, idate=idate, & dimx=lonb, dimy=latb, nfhour=nfhour, nfminute=nfminute, & nfsecondn=nfsecondn, nfsecondd=nfsecondd, iret=iret) else write(6,*)'***ERROR*** ','sfcf06mem001',' contains unrecognized format. ABORT' call MPI_Abort(MPI_COMM_WORLD,77,iret) stop endif endif nlon_ens = lonb nlat_ens = latb + 2 allocate(dtf_gsi(nlat_ens,nlon_ens),isli_gsi(nlat_ens,nlon_ens),work(nlat_ens,nlon_ens)) allocate(rwork1d(lonb*latb)) allocate(slmsk_ges(lonb,latb),slmsk_ens(lonb,latb)) allocate(dtf_ens(lonb,latb),dtzm(lonb,latb)) ! ! get slmsk for ges and anl of ensembles ! if ( sfcio ) then slmsk_ges = data_sfcges%slmsk ! read sfcgcy_mem001 and assign slmsk_ens call sfcio_srohdc(lun_sfcgcy,'sfcgcy_mem001',head_sfcgcy,data_sfcgcy,iret) slmsk_ens = data_sfcgcy%slmsk elseif ( nemsio ) then ! read land in fname_sfcges to get slmsk_ges call nemsio_readrecv(gfile_sfcges, 'land', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(mype,trim(my_name),'sfcf06_mem001','land','read',istop,iret) slmsk_ges=reshape(rwork1d,(/size(slmsk_ges,1),size(slmsk_ges,2)/)) ! close gfile_sfcges (not needed any more) call nemsio_close(gfile_sfcges, iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),'sfcf06_mem001',null,'close',istop,iret) ! open gfile_sfcgcy (sfcgcy_mem001) call nemsio_open(gfile_sfcgcy,'sfcgcy_mem001','read',iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),'sfcgcy_mem001',null,'open',istop,iret) ! read land in fname_sfcgcy to get slmsk_ens call nemsio_readrecv(gfile_sfcgcy, 'land', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(mype,trim(my_name),trim(fname_sfcgcy),'land','read',istop,iret) slmsk_ens=reshape(rwork1d,(/size(slmsk_ens,1),size(slmsk_ens,2)/)) endif ! get Tf analysis increment at ens resolution from anl resolution in GSI static analysis ! if ( (nlat_ens /= nlat_anl) .or. (nlon_ens /= nlon_anl) ) then if ( mype == 0 ) write(6,'(a,2(2a,2i5))')'getsfcnstensupdp: grid dimensions differ: ',& 'nlon_anl,nlat_anl = ',nlon_anl,nlat_anl,' nlon_ens,nlat_ens = ',nlon_ens,nlat_ens ! ! get lats and lons for GSI analysis grids ! jmax=nlat_anl-2 allocate(slatx(jmax),wlatx(jmax)) allocate(rlats_anl(nlat_anl),rlons_anl(nlon_anl)) call splat(idrt,jmax,slatx,wlatx) dlon=two*pi / real(nlon_anl,r_kind) do i=1,nlon_anl rlons_anl(i)=real(i-1,r_kind)*dlon enddo do i=1,(nlat_anl-1)/2 rlats_anl(i+1)=-asin(slatx(i)) rlats_anl(nlat_anl-i)=asin(slatx(i)) enddo rlats_anl(1)=-half*pi rlats_anl(nlat_anl)=half*pi deallocate(slatx,wlatx) ! get lats and lons for ensemble grids ! jmax=nlat_ens-2 allocate(slatx(jmax),wlatx(jmax)) allocate(rlats_ens(nlat_ens),rlons_ens(nlon_ens)) call splat(idrt,jmax,slatx,wlatx) dlon=two*pi / real(nlon_ens,r_kind) do i=1,nlon_ens rlons_ens(i)=real(i-1,r_kind)*dlon enddo do i=1,(nlat_ens-1)/2 rlats_ens(i+1)=-asin(slatx(i)) rlats_ens(nlat_ens-i)=asin(slatx(i)) enddo rlats_ens(1)=-half*pi rlats_ens(nlat_ens)=half*pi deallocate(slatx,wlatx) ! Get isli_gsi from slmsk_ens ! sumn = zero sums = zero do i=1, lonb sumn = slmsk_ens(i,1) + sumn sums = slmsk_ens(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 work(1,j)=sums do i=2, latb+1 work(i,j) = slmsk_ens(j,latb+2-i) end do work(latb+2,j)=sumn end do do j=1, nlon_ens do i=1, nlat_ens isli_gsi(i,j) = nint(work(i,j)) enddo enddo ! ! Get the expanded values for a surface type (0 = water now) and the new mask ! call int2_msk_glb_prep(dtf_anl,isli_anl,dtf_epd,isli_epd,nlat_anl,nlon_anl,0,nprep) ! ! Interpolate dtf_epd(nlat_anl,nlon_anl) to dtf_gsi(nlat_ens,nlon_ens) with surface mask accounted ! call int22_msk_glb(dtf_epd,isli_epd,rlats_anl,rlons_anl,nlat_anl,nlon_anl, & dtf_gsi,isli_gsi,rlats_ens,rlons_ens,nlat_ens,nlon_ens,0) ! ! transform the dtf_gsi(nlat_ens,nlon_ens) to dtf_ens(lonb,latb) for sfc file format ! do j=1,latb do i=1,lonb dtf_ens(i,j) = dtf_gsi(latb+2-j,i) enddo enddo else ! ! transform the dtf_gsi(nlat_ens,nlon_ens) to dtf_ens(lonb,latb) for sfc file format ! if ( mype == 0 ) write(6,'(a,2(a,2i5))')'getsfcnstensupdp: grid dimensions same: ',& 'nlon_anl,nlat_anl-2 = ',nlon_anl,nlat_anl-2,' lonb,latb = ',lonb,latb do j=1,latb do i=1,lonb dtf_ens(i,j)=dtf_anl(latb+2-j,i) enddo enddo endif ! if ( (nlat_ens /= nlat_anl) .or. (nlon_ens /= nlon_anl) ) then ! if ( npe < nanals ) then write(6,'(2(a,i5))')'***ERROR*** npe too small. npe = ',npe,' < nanals = ',nanals call MPI_Abort(MPI_COMM_WORLD,99,iret) stop endif mype1 = mype + 1 if ( mype1 > nanals ) then write (6,'(a,i5)') 'no files to process for mpi task = ',mype else ! get the filenames for each ensemble member write(charnanal,'(i3.3)') mype1 fname_nstges = 'nstf06_mem' // charnanal fname_nstanl = 'nstanl_mem' // charnanal fname_sfcgcy = 'sfcgcy_mem' // charnanal fname_sfcanl = 'sfcanl_mem' // charnanal ! ! update tsea & tref and then write out sfcanl & nstanl ! if (sfcio) then call sfcio_srohdc(lun_sfcgcy,trim(fname_sfcgcy),head_sfcgcy,data_sfcgcy,iret) if (mype==0) write(6,'(3a)')'Read ',trim(fname_sfcgcy),' in sfcio format ' call nstio_srohdc(lun_nstges,trim(fname_nstges),head_nst,data_nst,iret) if (mype==0) write(6,'(3a)')'Read ',trim(fname_nstges),' in sfcio format ' ! ! Assign sfcanl as sfcgcy ! head_sfcanl = head_sfcgcy data_sfcanl = data_sfcgcy ! ! For the new open water (sea ice just melted) grids, (1) set dtf_ens = zero (2) reset the NSSTM variables ! ! set tref = tfrozen = 271.2_r_kind ! note: data_sfcges%slmsk is the mask of the guess ! data_sfcanl%slmsk is the mask of the analysis ! where ( (slmsk_ens(:,:) == zero) .and. (slmsk_ges(:,:) == two) ) dtf_ens(:,:) = zero data_nst%xt(:,:) = zero data_nst%xs(:,:) = zero data_nst%xu(:,:) = zero data_nst%xv(:,:) = zero data_nst%xz(:,:) = z_w_max data_nst%zm(:,:) = zero data_nst%xtts(:,:) = zero data_nst%xzts(:,:) = zero data_nst%dt_cool(:,:) = zero data_nst%z_c(:,:) = zero data_nst%c_0(:,:) = zero data_nst%c_d(:,:) = zero data_nst%w_0(:,:) = zero data_nst%w_d(:,:) = zero data_nst%d_conv(:,:) = zero data_nst%ifd(:,:) = zero data_nst%tref(:,:) = tfrozen data_nst%qrain(:,:) = zero end where ! ! update analysis variable: Tref (foundation temperature) for nstanl file ! where ( slmsk_ens(:,:) == zero ) data_nst%tref(:,:) = max(data_nst%tref(:,:) + dtf_ens(:,:),tfrozen) elsewhere data_nst%tref(:,:) = data_sfcanl%tsea(:,:) end where ! Update guess date/time to analysis date/time for nst file head_nst%fhour = head_sfcanl%fhour ! forecast hour head_nst%idate(1) = head_sfcanl%idate(1) ! hour head_nst%idate(2) = head_sfcanl%idate(2) ! month head_nst%idate(3) = head_sfcanl%idate(3) ! day head_nst%idate(4) = head_sfcanl%idate(4) ! year ! Write updated information to nst analysis file call nstio_swohdc(lun_nstanl,trim(fname_nstanl),head_nst,data_nst,iret) write(6,101) trim(fname_nstanl),lonb,latb,head_nst%fhour,(head_nst%idate(i),i=1,4),iret 101 format(' sfcio_getsfcnstupdp: nst analysis written for ',& a,1x,2i6,1x,f4.1,4(i4,1x),' with iret = ',i5) ! ! update SST: tsea for sfcanl file ! if ( nst_gsi == 3 ) then call dtzm_2d(data_nst%xt,data_nst%xz,data_nst%dt_cool,data_nst%z_c, & data_sfcanl%slmsk,r_zsea1,r_zsea2,lonb,latb,dtzm) where ( slmsk_ens(:,:) == zero ) data_sfcanl%tsea(:,:) = max(data_nst%tref(:,:) + dtzm(:,:),tfrozen) end where ! Write updated information to surface analysis file call sfcio_swohdc(lun_sfcanl,trim(fname_sfcanl),head_sfcanl,data_sfcanl,iret) write(6,102) trim(fname_sfcanl),lonb,latb,head_sfcanl%fhour,(head_sfcanl%idate(i),i=1,4),iret 102 format(' sfcio_getsfcnstupdp: sfc analysis written for ',& a,1x,2i6,1x,f4.1,4(i4,1x),' with iret = ',i5) endif ! if ( nst_gsi == 3 ) then ! ! write out the info on the new open water and new sea ice grids ! if ( mype == 0 ) then n_new_water = 0 n_new_seaice = 0 do j=1,latb do i=1,lonb if ( slmsk_ens(i,j) == 0.0 .and. slmsk_ges(i,j) == 2.0 ) & n_new_water = n_new_water + 1 if ( slmsk_ens(i,j) == 2.0 .and. slmsk_ges(i,j) == 0.0 ) & n_new_seaice = n_new_seaice + 1 enddo enddo write(6,'(a,I3,2(1x,I8))') 'sfcio_getsfcnstens,nst_gsi,n_new_water,n_new_seaice:',nst_gsi,n_new_water,n_new_seaice endif elseif (nemsio) then if (mype==0) write(6,*)'computing mean with nemsio=',nemsio ! open gfile_sfcgcy (sfcgcy_mem???) call nemsio_open(gfile_sfcgcy,trim(fname_sfcgcy),'read',iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),'fname_sfcgcy',null,'open',istop,iret) ! open gfile_nstges (nstf06_mem???) call nemsio_open(gfile_nstges,trim(fname_nstges),'read',iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),'fname_nstges',null,'open',istop,iret) ! copy sfcgcy header info tosfcanl header, 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) if (iret /= 0) call error_msg(0,trim(my_name),trim(fname_sfcanl),null,'open',istop,iret) ! copy nstges header info to nstanl header gfile_nstanl=gfile_nstges ! open nemsio nstanl (with analysis time) call nemsio_open(gfile_nstanl,trim(fname_nstanl),'write',iret=iret,idate=idate, nfhour=nfhour,& nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd ) if (iret /= 0) call error_msg(0,trim(my_name),trim(fname_nstanl),null,'open',istop,iret) ! Allocate work array (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 end do ! ! read surface temperature from sfcgcy and save to tsea ! call nemsio_readrecv(gfile_sfcgcy, 'tmp', 'sfc', 1, rwork1d, iret=iret) if (iret /= 0) call error_msg(mype,trim(my_name),trim(fname_sfcgcy),'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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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(mype,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 where ( slmsk_ens(:,:) == zero .and. slmsk_ges(:,:) == two ) dtf_ens(:,:) = 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_ens(:,:) == zero ) tref(:,:) = max(tref(:,:) + dtf_ens(:,:),tfrozen) elsewhere tref(:,:) = tsea(:,:) end where ! ! update SST: tsea for sfc file with NSST profile ! if ( nst_gsi == 3 ) then 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_ens,r_zsea1,r_zsea2,lonb,latb,dtzm) where ( slmsk_ens(:,:) == zero ) tsea(:,:) = max(tref(:,:) + dtzm(:,:), tfrozen) end where 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(0,trim(my_name),trim(fname_sfcanl),'tmp','write',istop,iret) write(6,100) fname_sfcanl,lonb,latb,houra,idate(1:4),iret 100 format(' nemsio_getsfcnstensupdp: update tsea in sfcanl',a6,2i6,1x,f4.1,4(i4,1x),' with iret=',i2) ! ! update nsst records in nstanl ! ! slmsk rwork1d = reshape( slmsk_ens,(/size(rwork1d)/) ) call nemsio_writerecv(gfile_nstanl,'land','sfc',1,rwork1d,iret=iret) if (iret /= 0) call error_msg(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,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(0,trim(my_name),trim(fname_nstanl),'qrain','write',istop,iret) write(6,200) fname_nstanl,lonb,latb,houra,idate(1:4),iret 200 format(' nemsio_getsfcnstensupdp: update variables in nstanl',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_sfcgcy, iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),trim(fname_sfcgcy),null,'close',istop,iret) call nemsio_close(gfile_nstges, iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),trim(fname_nstges),null,'close',istop,iret) call nemsio_close(gfile_sfcanl,iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),trim(fname_sfcanl),null,'close',istop,iret) call nemsio_close(gfile_nstanl,iret=iret) if (iret /= 0) call error_msg(0,trim(my_name),trim(fname_nstanl),null,'close',istop,iret) endif ! if ( sfcio ) then endif ! if ( mype1 < nanals ) then call MPI_Barrier(MPI_COMM_WORLD,iret) if ( nproc == 0 ) call w3tage('GETSFCNSTENSUPDP') call MPI_Finalize(iret) if ( nproc == 0 .and. iret /= 0 ) & write(6,'(a,i5)'), 'MPI_Finalize error status, iret = ',iret contains subroutine error_msg(mype,sub_name,file_name,var_name,action,stop_code,error_code) use kinds, only: i_kind implicit none character(len=*), intent(in) :: sub_name,file_name,var_name,action integer(i_kind), intent(in) :: mype, stop_code, error_code if ( mype == 0 ) 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,'': problem opening file '',a,'', Status = '', i3)') & trim(sub_name), trim(file_name), error_code case('close') write(6,'(a,'': 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 ) then call MPI_Abort(MPI_COMM_WORLD,stop_code,iret) stop endif end subroutine error_msg END program getsfcnstensupdp