process_FVCOM.f90

Go to the documentation of this file.

program process_fvcom

   use mpi
   use kinds, only: r_kind, i_kind, r_single
   use module_ncio, only: ncio
   use module_nwp, only: fcst_nwp

   implicit none

! MPI variables
   integer :: npe, mype, mypelocal,ierror
   integer,allocatable :: mpi_layout_begin(:),mpi_layout_end(:)
!

!  New object-oriented declarations

   type(ncio) :: geo
   type(fcst_nwp) :: fcst

!  Grid variables

   character*180 :: thisfv3file
   character*2 :: workpath
   character*1 :: char1

   integer :: map_proj, nlon, nlat
   integer :: fv3lon, fv3lat, fv3times
   integer :: lbclon, lbclat, lbctimes
   integer :: i, j, t1, t2
   integer :: num_args, ix
   integer :: indexfvcomsel

   real :: rad2deg = 180.0/3.1415926
   real :: userdx, userdy, cen_lat, cen_lon
   real :: usertruelat1, usertruelat2, moad_cen_lat, stand_lon
   real :: truelat1, truelat2, stdlon, lat1, lon1, r_earth
   real :: knowni, knownj, dx
   real :: one, pi, deg2rad
   real :: zero

   character(len=180) :: fv3file
   character(len=180) :: fvcomfile
   character(len=180) :: wcstart
   character(len=180) :: inputfvcomselstr
   character(len=180), dimension(:), allocatable :: args
   integer            :: fv3_io_layout_y
   integer,allocatable :: fv3_layout_begin(:),fv3_layout_end(:)

   real(r_single), allocatable :: fv3ice(:,:), fv3sst(:,:)
   real(r_single), allocatable :: fv3sfct(:,:), fv3mask(:,:)
   real(r_single), allocatable :: fv3icet(:,:), fv3sfctl(:,:)
   real(r_single), allocatable :: fv3zorl(:,:), fv3hice(:,:)
   real(r_single), allocatable :: lbcice(:,:), lbcsst(:,:)
   real(r_single), allocatable :: lbcsfct(:,:), lbcmask(:,:)
   real(r_single), allocatable :: lbcicet(:,:), lbczorl(:,:)
   real(r_single), allocatable :: lbchice(:,:)

!  Declare namelists
!  SETUP (general control namelist) :

   integer :: update_type
   character(len=4) :: clayout_y
   integer :: iy,iblock

!   namelist/setup/update_type, t2

! MPI setup
   call mpi_init(ierror)
   call mpi_comm_size(mpi_comm_world,npe,ierror)
   call mpi_comm_rank(mpi_comm_world,mype,ierror)
!
! NCEP LSF has to use all cores allocated to run this application 
! but this if check can make sure only one core run through the real code.
!

   write(*,*) 'number of cores=',npe

   zero = 0.0
!  Get file names from command line arguements
   num_args = command_argument_count()
   allocate(args(num_args))

   do ix = 1, num_args
       call get_command_argument(ix,args(ix))
   end do
! fv3file: location of UFS grid
! fvcomfile: location of FVCOM data file
! wcstart: warm (restart) or cold start
! inputFVCOMtimes: string of time to use
   fv3file=trim(args(1))
   write(*,*) 'surface file=',trim(fv3file)
   fvcomfile=trim(args(2))    
   write(*,*) 'fvcom file=',trim(fvcomfile)
   wcstart=trim(args(3))
   write(*,*) 'warm or cold start = ', trim(wcstart)
   inputfvcomselstr=trim(args(4))
!   read(inputFVCOMselStr,*) inputFVCOMsel
   write(*,*) 'select time = ', trim(inputfvcomselstr)
   clayout_y=trim(args(5))
   read(clayout_y,'(I2)') fv3_io_layout_y
   if (wcstart == 'cold') then
     fv3_io_layout_y=1
   endif
   write(*,*) 'subdomain number=',fv3_io_layout_y
   if(mype < fv3_io_layout_y) then
      write(thisfv3file,'(a,I3.3)') 'stdout_fvcom.',mype
      open(6, file=trim(thisfv3file),form='formatted',status='unknown')
      write(6,*) '===> start process for core = ', mype
   endif
! 
! figure out subdomain grid begin and end grid index in Y for each core  
! The match the subdomains to each core
!
   allocate(fv3_layout_begin(fv3_io_layout_y))
   allocate(fv3_layout_end(fv3_io_layout_y))
   allocate(mpi_layout_begin(npe))
   allocate(mpi_layout_end(npe))
   fv3_layout_begin=0
   fv3_layout_end=0

   workpath='./'
   iy=0
   do ix=1,fv3_io_layout_y
      if(fv3_io_layout_y > 1) then
          write(thisfv3file,'(a,a,a1,I4.4)') trim(workpath), trim(fv3file),'.',ix-1
      else
          write(thisfv3file,'(a,a)') trim(workpath), trim(fv3file)
      endif
      call geo%open(trim(thisfv3file),'r',200)
      call geo%get_dim("yaxis_1",nlat)
      call geo%close
      fv3_layout_begin(ix)=iy+1
      fv3_layout_end(ix)=iy+nlat
      iy=fv3_layout_end(ix)
   enddo  ! find dimension
   write(6,'(a,20I5)') 'begin index for each subdomain',fv3_layout_begin
   write(6,'(a,20I5)') '  end index for each subdomain',fv3_layout_end

   mypelocal=mype+1
   if(npe==1) then
      mpi_layout_begin(mypelocal)=1
      mpi_layout_end(mypelocal)=fv3_io_layout_y
   else
      mpi_layout_begin=0
      mpi_layout_end=0
      do ix=1,fv3_io_layout_y
         iy=mod(ix-1,npe)+1
         mpi_layout_end(iy)=mpi_layout_end(iy)+1
      enddo
      iy=0
      do ix=1,npe
         if(mpi_layout_end(ix) > 0) then
           mpi_layout_begin(ix)=iy+1  
           mpi_layout_end(ix)=iy+mpi_layout_end(ix)
           iy=mpi_layout_end(ix)
         else
           mpi_layout_begin(ix)=0
           mpi_layout_end(ix)=0
         endif
      enddo 
   endif
   write(6,'(a)') 'begin and end domain index for each core:'
   write(6,'(20I5)') mpi_layout_begin
   write(6,'(20I5)') mpi_layout_end

if(mypelocal <= fv3_io_layout_y) then

   do ix=mpi_layout_begin(mypelocal),mpi_layout_end(mypelocal)

      write(6,*) 'process subdomain ',ix,' with core ',mype
      if(fv3_io_layout_y > 1) then
          write(thisfv3file,'(a,a,a1,I4.4)') trim(workpath), trim(fv3file),'.',ix-1
      else
          write(thisfv3file,'(a,a)') trim(workpath), trim(fv3file)
      endif
      write(6,*) 'sfc data file', trim(thisfv3file)

      call geo%open(trim(thisfv3file),'r',200)
      call geo%get_dim("xaxis_1",nlon)
      call geo%get_dim("yaxis_1",nlat)
      write(6,*) 'NLON,NLAT:', nlon, nlat
      call geo%close

      write(6,*) 'Finished reading sfc_data grid information.'
      write(6,*) ' '
!
!  figure out subdomain grid dimension in Y
!
      write(6,*) 'the Y dimensions=',fv3_layout_begin(ix),fv3_layout_end(ix)

!  Allocate variables for I/O

      allocate(fv3ice(nlon,nlat))
      allocate(fv3sfct(nlon,nlat))
      allocate(fv3sst(nlon,nlat))
      allocate(fv3mask(nlon,nlat))
      allocate(fv3icet(nlon,nlat))
      allocate(fv3sfctl(nlon,nlat))
      allocate(fv3zorl(nlon,nlat))
      allocate(fv3hice(nlon,nlat))

      allocate(lbcice(nlon,nlat))
      allocate(lbcsfct(nlon,nlat))
      allocate(lbcsst(nlon,nlat))
      allocate(lbcmask(nlon,nlat))
      allocate(lbcicet(nlon,nlat))
      allocate(lbczorl(nlon,nlat))
      allocate(lbchice(nlon,nlat))
!  Read fv3 sfc_data.nc before update

!   fv3file='sfc_data.nc'
!   fv3times: length of time dimension of UFS atmospheric grid (should be 1)
!   t1: index of time dimension to pull (should be 1)
      fv3times=1  
      t1=1

      call fcst%initial('FV3LAM',wcstart)
      call fcst%list_initial
      call fcst%read_n(trim(thisfv3file),'FV3LAM',wcstart,fv3lon,fv3lat,fv3times,t1,&
                   fv3mask,fv3sst,fv3ice,fv3sfct,fv3icet,fv3sfctl,fv3zorl,fv3hice,  &
                   1,nlat)
      call fcst%finish('FV3LAM',wcstart)


      write(6,*) 'fv3times: ', fv3times
      write(6,*) 'time to use: ', t1

!  Read FVCOM input datasets

!   fvcomfile='fvcom.nc'
!   lbctimes: length of time dimension of FVCOM input data (command line input)
! Space infront of ' FVCOM' below is important!!
      call fcst%initial(' FVCOM',wcstart)
      call fcst%list_initial
      call fcst%get_time_ind(trim(fvcomfile),inputfvcomselstr,indexfvcomsel)
!   t2: index of time dimension to pull from FVCOM
      t2=indexfvcomsel
      write(6,*) 'time asked for =', trim(inputfvcomselstr)
      write(6,*) 'time index selected = ', t2
      call fcst%read_n(trim(fvcomfile),' FVCOM',wcstart,lbclon,lbclat,lbctimes,t2, &
                   lbcmask,lbcsst,lbcice,lbcsfct,lbcicet,fv3sfctl,lbczorl,lbchice, &
                   fv3_layout_begin(ix),fv3_layout_end(ix))
      call fcst%finish(' FVCOM',wcstart)

!  Check that the dimensions match

      if (lbclon .ne. nlon .or. lbclat .ne. nlat) then
        write(6,*) 'ERROR: FVCOM/FV3 dimensions do not match:'
        write(6,*) 'lbclon: ', lbclon
        write(6,*) 'nlon: ', nlon
        write(6,*) 'lbclat: ', lbclat
        write(6,*) 'nlat: ', nlat
        stop 'error'
      endif

      write(6,*) 'lbctimes: ', lbctimes
      write(6,*) 'time to use: ', t2

      if (t2 .gt. lbctimes) then
        write(6,*) 'ERROR: Requested time dimension out of range'
        write(6,*) 'Length of time dimension: ', lbctimes
        write(6,*) 'Time index to use: ', t2
        stop 'error'
      endif

!  Update with FVCOM fields and process
!   ice cover data. Ice fraction is set
!   to a minimum of 15% due to FV3-LAM
!   raising any value below 15% to 15%.
   if (wcstart == 'warm') then
     do j=1,nlat
        do i=1,nlon
           if (lbcmask(i,j) > 0. .and. lbcsst(i,j) .ge. -90.0) then !GL Points
              !If ice fraction below 15%, set to 0
              if (lbcice(i,j) < 0.15) then
                lbcice(i,j) = 0.0
                lbchice(i,j) = 0.0 !remove ice thickness
              endif
              fv3ice(i,j) = lbcice(i,j)
              fv3hice(i,j) = lbchice(i,j)

              !If ice in FVCOM, but not in FV3-LAM, change to ice
              if (lbcice(i,j) > 0. .and. fv3mask(i,j) == 0.) then
                fv3mask(i,j) = 2.
                fv3zorl(i,j) = 1.1
              endif
              !If ice in FV3-LAM and not FVCOM, remove it from FV3-LAM
              if (fv3mask(i,j) == 2. .and. lbcice(i,j) == 0.) then
                fv3mask(i,j) = 0.
                fv3zorl(i,j) = zero / zero !Use Fill_Value
              endif
              fv3sst(i,j) = lbcsst(i,j) + 273.15
              fv3sfct(i,j) = lbcsst(i,j) + 273.15
              fv3icet(i,j) = lbcsst(i,j) + 273.15
              fv3sfctl(i,j)= lbcsst(i,j) + 273.15
               !If ice exists in FVCOM, change ice surface temp
              if (lbcice(i,j) > 0.) then
                fv3icet(i,j) = lbcicet(i,j) + 273.15
              end if
           end if
        enddo
     enddo
   else if (wcstart == 'cold') then
     do j=1,nlat
        do i=1,nlon
           if (lbcmask(i,j) > 0. .and. lbcsst(i,j) .ge. -90.0) then
              !If ice fraction below 15%, set to 0
              if (lbcice(i,j) < 0.15) then
                lbcice(i,j) = 0.0
                lbchice(i,j) = 0.0 !remove ice thickness
              endif
              fv3ice(i,j) = lbcice(i,j)
              fv3hice(i,j) = lbchice(i,j)
              !If ice in FVCOM, but not in FV3-LAM, change to ice
              if (lbcice(i,j) > 0. .and. fv3mask(i,j) == 0.) then
                fv3mask(i,j) = 2.
                fv3zorl(i,j) = 1.1
              endif
              !If ice in FV3-LAM and not FVCOM, remove it from FV3-LAM
              if (fv3mask(i,j) == 2. .and. lbcice(i,j) == 0.) then
                fv3mask(i,j) = 0.
                fv3zorl(i,j) = zero / zero !Use Fill_Value
              endif
              fv3sst(i,j) = lbcsst(i,j) + 273.15
              fv3sfct(i,j) = lbcsst(i,j) + 273.15
              fv3icet(i,j) = lbcsst(i,j) + 273.15
               !If ice exists in FVCOM, change ice surface temp
              if (lbcice(i,j) > 0.) then
                fv3icet(i,j) = lbcicet(i,j) + 273.15
              end if
           end if
        enddo
     enddo
   else
     write(6,*) 'Variable wcstart is not set to either warm or cold'
   end if

! Write out sfc file again

      write(6,*) "udpate file=",trim(thisfv3file)
      call geo%open(trim(thisfv3file),'w',300)
      call geo%replace_var("tsea",nlon,nlat,fv3sst)
      call geo%replace_var("fice",nlon,nlat,fv3ice)
      call geo%replace_var("slmsk",nlon,nlat,fv3mask)
      call geo%replace_var("tisfc",nlon,nlat,fv3icet)
      call geo%replace_var("hice",nlon,nlat,fv3hice)
      if (wcstart == 'cold') then
! Add_New_Var takes names of (Variable,Dim1,Dim2,Dim3,Long_Name,Units)
         call geo%replace_var("zorl",nlon,nlat,fv3zorl)
         call geo%add_new_var('glmsk','xaxis_1','yaxis_1','Time','glmsk','none','float')
         call geo%replace_var('glmsk',nlon,nlat,lbcmask)
      end if
      if (wcstart == 'warm') then
         call geo%replace_var("zorli",nlon,nlat,fv3zorl)
         call geo%replace_var("tsfc",nlon,nlat,fv3sfct)
         call geo%replace_var("tsfcl",nlon,nlat,fv3sfctl)
         call geo%add_new_var('glmsk','xaxis_1','yaxis_1','Time','glmsk','none','float')
         call geo%replace_var('glmsk',nlon,nlat,lbcmask)
      end if
      call geo%close

      deallocate(fv3ice)
      deallocate(fv3sfct)
      deallocate(fv3sst)
      deallocate(fv3mask)
      deallocate(fv3icet)
      deallocate(fv3sfctl)
      deallocate(fv3zorl)
      deallocate(fv3hice)

      deallocate(lbcice)
      deallocate(lbcsfct)
      deallocate(lbcsst)
      deallocate(lbcmask)
      deallocate(lbcicet)
      deallocate(lbczorl)
      deallocate(lbchice)

      write(6,*) "=== LOWBC UPDATE SUCCESS ===",ix

   enddo ! loop through subdomain

endif ! mype==0

call mpi_finalize(ierror)


end program process_fvcom