program nstrtg
!
!$$$  main program documentation block
! Abstract:  Generate a RTG-like SST analysis file at a lat/lon grids
!           (resolution defined by namelist variables)
!           
! Notes of the procudure to generate RTG-like SST file from NSST Tf analysis :
!       1. Read namelist variables
!       2. Get lat/lon for the target grids (ny) (nx)
!       3. Read surface mask at the target grids (nx,ny)
!       4. Get 6-hourly GFS NSST foundation temperature analysis (T1534 Gaussian grids, at prsent))
!       5. Get sea ice analysis at target grids (interpolated when necessary)
!       6. Modify the mask with sea ice info
!       7. get tf_clm at the analysis time & at target resolution from monthly climatology (nx,ny)
!       8. get sal_clm at the analysis time & at target resolution from monthly climatology (nx,ny)
!       9. interpolate tf_nst (Gaussian lat/lon) to tf (lat/lon)
!      10. set tf to be SST climatology for the missing water grids
!      11. get water temperature for grids with positive sea ice fraction 
!      12. fill the land grids temperature with known water grids temperatures
!      13. Write tf as grib format w/o bitmap
!     
! Created by Xu Li, Mar., 2019

!$$$

 use set_para, only: init_setup

 implicit none

 integer, parameter :: sfcflag=0,miss_fill=0
 real,    parameter :: bmiss=-999.0, tzero=273.16,tfrozen=271.20
!
! Arrays at Gaussian grids
!
 real,    allocatable, dimension(:,:) :: tf_nst              ! NSST Tf analysis from nemsio file (Gaussian grids)
 integer, allocatable, dimension(:,:) :: mask_nst            ! mask of NSST Tf (0 = water, others = non-water)
 real,    allocatable, dimension(:)   :: xlats_nst           ! latitudes  of NSST Tf at Gaussian grids
 real,    allocatable, dimension(:)   :: xlons_nst           ! longitudes of NSST Tf at Gaussian grids
!
! Arrays at target grids (nx,ny), (nx) or (ny)
!
 real,    allocatable, dimension(:,:) :: tf                  ! NSST Tf at target grids (nx,ny)
 real,    allocatable, dimension(:,:) :: sice                ! Sea ice analysis at target grids (nx,ny)
 real,    allocatable, dimension(:,:) :: sal_clm             ! Salinity Climatology at target grids (nx,ny)
 real,    allocatable, dimension(:,:) :: tf_clm              ! SST Climatology at target grids and valid time (nx,ny)
 real,    allocatable, dimension(:)   :: xlats               ! latitudes of target grids (ny)
 real,    allocatable, dimension(:)   :: xlons               ! longitudes of target grids (nx)
 integer, allocatable, dimension(:,:) :: mask                ! mask at target grids (0 = water, 1 = land, 2 = sea ice) (nx,ny)
 integer, allocatable, dimension(:,:) :: masks               ! mask at target grids (0 = water or sea ice , 1 = land) (nx,ny)
!
! Arrays target grids (nx*ny)
!
 real,    allocatable, dimension(:)   :: tf_ij               ! NSST Tf at target grids (nx*ny)
 real,    allocatable, dimension(:)   :: xlats_ij            ! latitudes of target grids (nx*ny)
 real,    allocatable, dimension(:)   :: xlons_ij            ! longitudes of target grids (nx*ny)
 integer, allocatable, dimension(:)   :: mask_ij             ! mask at target grids (0 =water, 1 = land, 2 = sea ice) (nx*ny)
 integer, allocatable, dimension(:)   :: masks_ij            ! mask at target grids (0 = water or sea ice , 1 = land) (nx*ny)

 real :: tfreez

!input/output data file names
 character (len=6),  parameter :: fin_tf_clm='sstclm'        ! SST climatology 
 character (len=6),  parameter :: fin_sfcanl='sfcanl'        ! GFS surface analysis file, 6-hourly
 character (len=6),  parameter :: fin_iceanl='iceanl'        ! NCEP se ice analysis file, daily
 character (len=6),  parameter :: fin_salclm='salclm'        ! Levitus Salinity Climatology (one degree)
 character (len=6),  parameter :: fin_rtgmsk='rtgmsk'        ! land/water mask

 character (len=10), parameter :: fout_tf_grb='tf_rtg_grb'
 character (len=16), parameter :: fout_tf_grb_awips='tf_rtg_grb_awips'


 integer :: lensfc,mon1,mon2
 integer :: iy,im,id,ih,i,j,k,ij
 integer :: nxc,nyc,nxc_sal,nyc_sal
 integer :: nx_nst,ny_nst,nx,ny
 integer :: nfill_clm_water,nfill_clm_land,nfill_clm_ice
 real    :: wei1,wei2,dsearch,x0,y0,dres
 integer :: maskflag      ! 1 = with bitmap; 0 = w/o bitmap
 logical :: lputsi
 character (len=10) :: catime
 namelist/setup/catime,lputsi,dsearch,nx,ny

 call init_setup

 CALL W3TAGB('nsst_tf_to_rtg',2019,0,0,'emc')

!Read namelist

 read(5,setup)

 read(catime(1:10),'(i4,3i2)') iy,im,id,ih

 allocate( tf(nx,ny),sice(nx,ny),tf_clm(nx,ny),sal_clm(nx,ny),mask(nx,ny),masks(nx,ny) )
 allocate( mask_ij(nx*ny),masks_ij(nx*ny),tf_ij(nx*ny))
 allocate( xlats(ny),xlons(nx),xlats_ij(nx*ny),xlons_ij(nx*ny) )
!
! get xlats & xlons
!
  dres = 180.0/real(ny)
  y0 = 0.5*dres-90.0
  x0 = 0.5*dres

! Get lat_sst & lon_sst
  do j = 1, ny
     xlats(j) = y0 + real(j-1)*dres
  enddo

  do i = 1, nx
     xlons(i) = (x0 + real(i-1)*dres)
  enddo
!
! Get 1-dimentional xlats and xlons for global grids (nx*ny)
!
 ij = 0
 do j = 1, ny
    do i = 1, nx
       ij = ij + 1
       xlats_ij(ij) = xlats(j)
       xlons_ij(ij) = xlons(i)
    enddo
 enddo
!
! read target mask (nx,ny): mask = 0 for ocean; mask = 1 for land
!
 call get_mask(fin_rtgmsk,mask,nx,ny)
 masks = mask
 write(*,*) 'get mask from ',fin_rtgmsk
!
! Get dimensions of sfcanl, allocate arrays and read NSST tf
!
 call get_sfcanl_dim(fin_sfcanl,ny_nst,nx_nst)
 write(*,'(a,8I5)') 'iy,im,id,ih,ny_nst,nx_nst,ny,nx : ',iy,im,id,ih,ny_nst,nx_nst,ny,nx
 allocate( tf_nst(nx_nst,ny_nst),mask_nst(nx_nst,ny_nst) )
 allocate( xlons_nst(nx_nst),xlats_nst(ny_nst) )
 call read_tf_sfcanl(fin_sfcanl,tf_nst,mask_nst,xlats_nst,xlons_nst,ny_nst,nx_nst)
!
! get sea ice analysis (daily and available 6-houly files) at taget grids (nx,ny)
!
  call get_seaice(xlats_ij,xlons_ij,ny,nx,sice)
!
! modify mask for sea ice
!
 do j = 1, ny
    do i = 1, nx
       if ( sice(i,j) > 0.50 ) then
          mask(i,j) = 2
          masks(i,j) = 0
       endif
    enddo
 enddo
 lensfc=nx*ny  
 mask_ij(:)  = reshape( mask,  (/lensfc/) )
 masks_ij(:) = reshape( masks, (/lensfc/) )
!
! get tf_clm at the analysis time & at target resolution from monthly climatology (nx,ny)
!
  call get_tf_clm(mask_ij,xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,dsearch,miss_fill,bmiss,tf_clm)
!
! get sal_clm at the analysis time & at target resolution from monthly climatology (nx,ny)
!
 call get_sal_clm(mask_ij,xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,dsearch,miss_fill,bmiss,sal_clm)
!
! interpolate tf_nst (Gaussian lat/lon) to tf (lat/lon)
!
 call lalo_to_tile(tf_nst, mask_nst,xlats_nst,xlons_nst,ny_nst,nx_nst, &
                   tf_ij,  mask_ij, xlats_ij, xlons_ij, ny,nx, &
                   sfcflag,dsearch, miss_fill,bmiss)
 write(*,*) 'done with lalo_to_tile for tf_ij'

 tf(:,:) = reshape (tf_ij, (/nx,ny/) )
!
! set temp to be RTG SST climatology for the missing water grids 
!
 nfill_clm_water=0
 nfill_clm_land=0
 nfill_clm_ice=0
 do j = 1, ny
    do i = 1, nx
       if ( (mask(i,j) == 0 .and. tf(i,j) == bmiss) .or. mask(i,j) /= 0 ) then
          tf(i,j) = tf_clm(i,j)
!
!         get the counts statistics for a surface type
!
          if ( mask(i,j) == 0 ) then
             nfill_clm_water = nfill_clm_water + 1
          elseif ( mask(i,j) == 1 ) then
             nfill_clm_land = nfill_clm_land + 1
          else
             nfill_clm_ice = nfill_clm_ice + 1
          endif
       endif
    enddo
 enddo

 write(*,*) 'nfill_clm_water : ',nfill_clm_water
 write(*,*) 'nfill_clm_land : ',nfill_clm_land
 write(*,*) 'nfill_clm_ice : ',nfill_clm_ice
!
! calculate water temperature for grids with positive sea ice fraction
! set temperature of some fresh water lakes to 0C
!
 do j = 1, ny
    do i = 1, nx
       if ( mask(i,j) == 2 ) then
!         great lakes
          if ( xlats(j) > 41.0  .and. xlats(j) < 50.75 .and. &
               xlons(i) > 266.0 .and. xlons(i) < 285.0 ) then
             tf(i,j) = tzero
!         lake winnepeg, great bear lake, great slave lake, et al
          elseif ( xlats(j) > 50.0  .and. xlats(j) < 67.0 .and. &
                   xlons(i) > 235.0 .and. xlons(i) < 265.0 ) then
              tf(i,j) = tzero
          else
             if ( sal_clm(i,j) < 5.0 .or. sal_clm(i,j) > 50.0 ) then     ! non-water
                tf(i,j) = tfreez(sal_clm(i,j))+tzero
             endif
          endif
       endif
    enddo
 enddo

!
! fill land grids temperature with known water grids temperatures
!
 if ( lputsi ) then
    tf(:,:) = tf(:,:) - tzero
    call putsi(tf,mask,nx,ny)
    tf(:,:) = tf(:,:) + tzero
 endif
!
! write out tf as grib1 and grib2 format
!
 maskflag=0      ! no bitmap
 call togrib(fout_tf_grb,tf,mask,maskflag,nx,ny,iy,im,id,ih)

 maskflag=1      ! with bitmap
 call togrib(fout_tf_grb_awips,tf,mask,maskflag,nx,ny,iy,im,id,ih)

 call w3tage('nsst_tf_to_rtg')

 write(*,*) 'All Done'

end program nstrtg