cycle.f90

Go to the documentation of this file.

  
 PROGRAM sfc_drv

 use mpi

 IMPLICIT NONE
!
 CHARACTER(LEN=3) :: donst
 INTEGER :: idim, jdim, lsoil, lugb, iy, im, id, ih, ialb
 INTEGER :: isot, ivegsrc, lensfc, zsea1_mm, zsea2_mm, ierr
 INTEGER :: nprocs, myrank, num_threads, num_parthds, max_tasks
 REAL    :: fh, deltsfc, zsea1, zsea2
 LOGICAL :: use_ufo, do_nsst, do_lndinc, do_sfccycle, frac_grid
!
 NAMELIST/namcyc/ idim,jdim,lsoil,lugb,iy,im,id,ih,fh,&
                  deltsfc,ialb,use_ufo,donst,             &
                  do_sfccycle,isot,ivegsrc,zsea1_mm,      &
                  zsea2_mm, max_tasks, do_lndinc, frac_grid
!
 DATA idim,jdim,lsoil/96,96,4/
 DATA iy,im,id,ih,fh/1997,8,2,0,0./
 DATA lugb/51/, deltsfc/0.0/, ialb/1/, max_tasks/99999/
 DATA isot/1/, ivegsrc/2/, zsea1_mm/0/, zsea2_mm/0/
!
 CALL mpi_init(ierr)
 CALL mpi_comm_size(mpi_comm_world, nprocs, ierr)
 CALL mpi_comm_rank(mpi_comm_world, myrank, ierr)

 if (myrank==0) call w3tagb('GLOBAL_CYCLE',2018,0179,0055,'NP20')

 num_threads = num_parthds()

 print*
 print*,"STARTING CYCLE PROGRAM ON RANK ", myrank
 print*,"RUNNING WITH ", nprocs, "TASKS"
 print*,"AND WITH ", num_threads, " THREADS."

 use_ufo = .false.
 donst   = "NO"
 do_lndinc   = .false.
 do_sfccycle = .true.
 frac_grid = .false.

 print*
 print*,"READ NAMCYC NAMELIST."

 CALL baopenr(36, "fort.36", ierr)
 READ(36, nml=namcyc)
!IF (MYRANK==0) WRITE(6,NAMCYC)

 IF (max_tasks < 99999 .AND. myrank > (max_tasks - 1)) THEN
   print*,"USER SPECIFIED MAX NUMBER OF TASKS: ", max_tasks
   print*,"WILL NOT RUN CYCLE PROGRAM ON RANK: ", myrank
   GOTO 333
 ENDIF

 lensfc = idim*jdim ! TOTAL NUMBER OF POINTS FOR THE CUBED-SPHERE TILE

 zsea1 = float(zsea1_mm) / 1000.0  ! CONVERT FROM MM TO METERS
 zsea2 = float(zsea2_mm) / 1000.0

 IF (donst == "YES") THEN
   do_nsst=.true.
 ELSE
   do_nsst=.false.
 ENDIF

 print*
 IF (myrank==0) print*,"LUGB,IDIM,JDIM,ISOT,IVEGSRC,LSOIL,DELTSFC,IY,IM,ID,IH,FH: ", &
              lugb,idim,jdim,isot,ivegsrc,lsoil,deltsfc,iy,im,id,ih,fh

 CALL sfcdrv(lugb,idim,jdim,lensfc,lsoil,deltsfc,  &
             iy,im,id,ih,fh,ialb,                  &
             use_ufo,do_nsst,do_sfccycle,do_lndinc, &
             frac_grid,zsea1,zsea2,isot,ivegsrc,myrank)
 
 print*
 print*,'CYCLE PROGRAM COMPLETED NORMALLY ON RANK: ', myrank

 333 CONTINUE

 CALL mpi_barrier(mpi_comm_world, ierr)

 if (myrank==0) call w3tage('GLOBAL_CYCLE')

 CALL mpi_finalize(ierr)

 stop

 END PROGRAM sfc_drv

 SUBROUTINE sfcdrv(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
                   IY,IM,ID,IH,FH,IALB,                  &
                   USE_UFO,DO_NSST,DO_SFCCYCLE,DO_LNDINC,&
                   FRAC_GRID,ZSEA1,ZSEA2,ISOT,IVEGSRC,MYRANK)
!
 USE read_write_data
 use machine
 USE mpi
 USE land_increments, ONLY: add_increment_soil,     &
                            add_increment_snow,     &
                            calculate_landinc_mask, &
                            apply_land_da_adjustments_soil, &
                            apply_land_da_adjustments_snd, &
                            lsm_noah, lsm_noahmp

 IMPLICIT NONE

 INTEGER, INTENT(IN) :: IDIM, JDIM, LENSFC, LSOIL, IALB
 INTEGER, INTENT(IN) :: LUGB, IY, IM, ID, IH
 INTEGER, INTENT(IN) :: ISOT, IVEGSRC, MYRANK

 LOGICAL, INTENT(IN) :: USE_UFO, DO_NSST,DO_SFCCYCLE
 LOGICAL, INTENT(IN) :: DO_LNDINC, FRAC_GRID
 
 REAL, INTENT(IN)    :: FH, DELTSFC, ZSEA1, ZSEA2

 INTEGER, PARAMETER  :: NLUNIT=35
 INTEGER, PARAMETER  :: SZ_NML=1

 CHARACTER(LEN=5)    :: TILE_NUM
 CHARACTER(LEN=500)  :: NST_FILE
 CHARACTER(LEN=500)  :: LND_SOI_FILE
 CHARACTER(LEN=4)    :: INPUT_NML_FILE(SZ_NML)

 INTEGER             :: I, IERR
 INTEGER             :: I_INDEX(LENSFC), J_INDEX(LENSFC)
 INTEGER             :: IDUM(IDIM,JDIM)
 integer             :: num_parthds, num_threads

 LOGICAL             :: IS_NOAHMP
 INTEGER             :: LSM

 real(kind=kind_io8) :: min_ice(lensfc)

 REAL                :: SLMASK(LENSFC), OROG(LENSFC)
 REAL                :: SIHFCS(LENSFC), SICFCS(LENSFC)
 REAL                :: SITFCS(LENSFC), TSFFCS(LENSFC)
 REAL                :: SWEFCS(LENSFC), ZORFCS(LENSFC)
 REAL                :: ALBFCS(LENSFC,4), TG3FCS(LENSFC)
 REAL                :: CNPFCS(LENSFC), SMCFCS(LENSFC,LSOIL)
 REAL                :: STCFCS(LENSFC,LSOIL), SLIFCS(LENSFC)
 REAL                :: AISFCS(LENSFC), F10M(LENSFC)
 REAL                :: VEGFCS(LENSFC), VETFCS(LENSFC)
 REAL                :: SOTFCS(LENSFC), ALFFCS(LENSFC,2)
 REAL                :: CVFCS(LENSFC), CVTFCS(LENSFC)
 REAL                :: CVBFCS(LENSFC), TPRCP(LENSFC)
 REAL                :: SRFLAG(LENSFC), SNDFCS(LENSFC)
 REAL                :: SLCFCS(LENSFC,LSOIL), VMXFCS(LENSFC)
 REAL                :: VMNFCS(LENSFC), T2M(LENSFC)
 REAL                :: Q2M(LENSFC), SLPFCS(LENSFC)
 REAL                :: ABSFCS(LENSFC), OROG_UF(LENSFC)
 REAL                :: USTAR(LENSFC)
 REAL                :: FMM(LENSFC), FHH(LENSFC)
 REAL                :: RLA(LENSFC), RLO(LENSFC)
 REAL(KIND=4)        :: ZSOIL(LSOIL)
 REAL                :: SIG1T(LENSFC)
 REAL, ALLOCATABLE   :: STC_BCK(:,:), SMC_BCK(:,:), SLC_BCK(:,:)
 REAL, ALLOCATABLE   :: SLIFCS_FG(:), SICFCS_FG(:)
 INTEGER, ALLOCATABLE :: LANDINC_MASK_FG(:), LANDINC_MASK(:)
 REAL, ALLOCATABLE   :: SND_BCK(:), SND_INC(:), SWE_BCK(:)
 REAL(KIND=KIND_IO8), ALLOCATABLE :: SLMASKL(:), SLMASKW(:), LANDFRAC(:)

 TYPE(NSST_DATA)     :: NSST
 real, dimension(idim,jdim) :: tf_clm,tf_trd,sal_clm
 real, dimension(lensfc)    :: tf_clm_tile,tf_trd_tile,sal_clm_tile
 INTEGER             :: veg_type_landice
 INTEGER, DIMENSION(LENSFC) :: STC_UPDATED, SLC_UPDATED

 LOGICAL :: FILE_EXISTS, DO_SOI_INC, DO_SNO_INC
 CHARACTER(LEN=3)       :: RANKCH

!--------------------------------------------------------------------------------
! NST_FILE is the path/name of the gaussian GSI file which contains NSST
! increments.
!--------------------------------------------------------------------------------
 
 NAMELIST/namsfcd/ nst_file, lnd_soi_file, do_sno_inc

 DATA nst_file/'NULL'/
 DATA lnd_soi_file/'NULL'/

 do_sno_inc = .false.
 do_soi_inc = .false.
 
 sig1t = 0.0            ! Not a dead start!

 input_nml_file = "NULL"

 CALL baopenr(37, "fort.37", ierr)
 READ (37, nml=namsfcd)

 print*
 print*,'IN ROUTINE SFCDRV,IDIM=',idim,'JDIM=',jdim,'FH=',fh

!--------------------------------------------------------------------------------
! READ THE OROGRAPHY AND GRID POINT LAT/LONS FOR THE CUBED-SPHERE TILE.
!--------------------------------------------------------------------------------

 ALLOCATE(landfrac(lensfc))
 IF(frac_grid) THEN
   print*,'- RUNNING WITH FRACTIONAL GRID.'
   CALL read_lat_lon_orog(rla,rlo,orog,orog_uf,tile_num,idim,jdim,lensfc,landfrac=landfrac)
 ELSE
   CALL read_lat_lon_orog(rla,rlo,orog,orog_uf,tile_num,idim,jdim,lensfc)
   landfrac=-999.9
 ENDIF

 DO i = 1, idim
   idum(i,:) = i
 ENDDO

 i_index = reshape(idum, (/lensfc/))

 DO i = 1, jdim
   idum(:,i) = i
 ENDDO

 j_index = reshape(idum, (/lensfc/) )

 IF (do_nsst) THEN
   print*
   print*,"WILL PROCESS NSST RECORDS."
   ALLOCATE(nsst%C_0(lensfc))
   ALLOCATE(nsst%C_D(lensfc))
   ALLOCATE(nsst%D_CONV(lensfc))
   ALLOCATE(nsst%DT_COOL(lensfc))
   ALLOCATE(nsst%IFD(lensfc))
   ALLOCATE(nsst%QRAIN(lensfc))
   ALLOCATE(nsst%TREF(lensfc))
   ALLOCATE(nsst%TFINC(lensfc))
   ALLOCATE(nsst%W_0(lensfc))
   ALLOCATE(nsst%W_D(lensfc))
   ALLOCATE(nsst%XS(lensfc))
   ALLOCATE(nsst%XT(lensfc))
   ALLOCATE(nsst%XTTS(lensfc))
   ALLOCATE(nsst%XU(lensfc))
   ALLOCATE(nsst%XV(lensfc))
   ALLOCATE(nsst%XZ(lensfc))
   ALLOCATE(nsst%XZTS(lensfc))
   ALLOCATE(nsst%Z_C(lensfc))
   ALLOCATE(nsst%ZM(lensfc))
   ALLOCATE(slifcs_fg(lensfc))
   ALLOCATE(sicfcs_fg(lensfc))
 ENDIF

IF (do_lndinc) THEN
   IF (frac_grid) THEN
     print *, 'FATAL ERROR: land increment update does not work with fractional grids.'
     call mpi_abort(mpi_comm_world, 17, ierr)
   ENDIF
   ! identify variables to be updated, and allocate arrays.
   IF  (trim(lnd_soi_file) .NE. "NULL") THEN
       do_soi_inc = .true.
       print*
       print*," APPLYING SOIL INCREMENTS FROM THE GSI"
       ALLOCATE(stc_bck(lensfc, lsoil), smc_bck(lensfc, lsoil), slc_bck(lensfc,lsoil))
       ALLOCATE(landinc_mask_fg(lensfc))
   ENDIF
   ! FOR NOW, CODE SO CAN DO BOTH, BUT MIGHT NEED TO THINK ABOUT THIS SOME MORE.
   IF  (do_sno_inc) THEN
       ! ideally, would check here that sfcsub snow DA update is not also requested
       ! but latter is controlled by fnsol, which is read in within that routine.
       ! should be done at script level.
       print*
       print*," APPLYING SNOW INCREMENTS FROM JEDI"
       ALLOCATE(snd_bck(lensfc), snd_inc(lensfc), swe_bck(lensfc))
   ENDIF
   ! set-up land mask info
   ALLOCATE(landinc_mask(lensfc))
   if (ivegsrc == 2) then   ! sib
      veg_type_landice=13
   else
      veg_type_landice=15
   endif
ENDIF

!--------------------------------------------------------------------------------
! READ THE INPUT SURFACE DATA ON THE CUBED-SPHERE TILE.
!--------------------------------------------------------------------------------

 CALL read_data(lsoil,lensfc,do_nsst,.false.,is_noahmp=is_noahmp, &
                tsffcs=tsffcs,smcfcs=smcfcs,   &
                swefcs=swefcs,stcfcs=stcfcs,tg3fcs=tg3fcs,zorfcs=zorfcs,  &
                cvfcs=cvfcs,  cvbfcs=cvbfcs,cvtfcs=cvtfcs,albfcs=albfcs,  &
                vegfcs=vegfcs,slifcs=slifcs,cnpfcs=cnpfcs,f10m=f10m    ,  &
                vetfcs=vetfcs,sotfcs=sotfcs,alffcs=alffcs,ustar=ustar  ,  &
                fmm=fmm      ,fhh=fhh      ,sihfcs=sihfcs,sicfcs=sicfcs,  &
                sitfcs=sitfcs,tprcp=tprcp  ,srflag=srflag,sndfcs=sndfcs,  &
                vmnfcs=vmnfcs,vmxfcs=vmxfcs,slcfcs=slcfcs,slpfcs=slpfcs,  &
                absfcs=absfcs,t2m=t2m      ,q2m=q2m      ,slmask=slmask,  &
                zsoil=zsoil,   nsst=nsst)

 IF (frac_grid .AND. .NOT. is_noahmp) THEN
   print *, 'FATAL ERROR: NOAH lsm update does not work with fractional grids.'
   call mpi_abort(mpi_comm_world, 18, ierr)
 ENDIF

 IF (frac_grid .AND. do_sno_inc) THEN
   print *, 'FATAL ERROR: Snow increment update does not work with fractional grids.'
   call mpi_abort(mpi_comm_world, 19, ierr)
 ENDIF

 IF (is_noahmp .AND. do_sno_inc) THEN
   print *, 'FATAL ERROR: Snow increment update does not work with NOAH_MP.'
   call mpi_abort(mpi_comm_world, 29, ierr)
 ENDIF

 IF (is_noahmp) THEN 
        lsm=lsm_noahmp
 ELSE
        lsm=lsm_noah
 ENDIF

 IF (use_ufo) THEN
   print*
   print*,'USE UNFILTERED OROGRAPHY.'
 ELSE
   orog_uf = 0.0
 ENDIF
 
 DO i=1,lensfc
   aisfcs(i) = 0.
   IF(nint(slifcs(i)).EQ.2) aisfcs(i) = 1.
 ENDDO

 IF (do_nsst) THEN
   sicfcs_fg=sicfcs
   IF (.NOT. do_sfccycle ) THEN
     print*
     print*,"FIRST GUESS MASK ADJUSTED BY IFD RECORD"
     slifcs_fg = slifcs
     WHERE(nint(nsst%IFD) == 3) slifcs_fg = 2.0
   ELSE
     print*
     print*,"SAVE FIRST GUESS MASK"
     slifcs_fg = slifcs
   ENDIF
 ENDIF
 
 ! CALCULATE MASK FOR LAND INCREMENTS
 IF (do_lndinc)  &
    CALL calculate_landinc_mask(smcfcs(:,1),swefcs, vetfcs,  &
                    lensfc,veg_type_landice,  landinc_mask)

!--------------------------------------------------------------------------------
! UPDATE SURFACE FIELDS.
!
! FIRST, SET WATER AND LAND MASKS - SLMASKW/SLMASKL. FOR UNCOUPLED
! (NON-FRACTIONAL) MODE, THESE ARE IDENTICAL TO THE CURRENT
! MASK - '0' WATER; '1' LAND.
!--------------------------------------------------------------------------------

 IF (do_sfccycle) THEN
   ALLOCATE(slmaskl(lensfc), slmaskw(lensfc))

   set_mask : IF (frac_grid) THEN

     DO i=1,lensfc
       IF(landfrac(i) > 0.0_kind_io8) THEN
         slmaskl(i) = ceiling(landfrac(i)-1.0e-6_kind_io8)
         slmaskw(i) =   floor(landfrac(i)+1.0e-6_kind_io8)
       ELSE
         IF(nint(slmask(i)) == 1) THEN ! If landfrac is zero, this should not happen.
                                       ! So, stop processing.
           print*, 'FATAL ERROR: LAND FRAC AND SLMASK MISMATCH.'
           CALL mpi_abort(mpi_comm_world, 27, ierr)
         ELSE
           slmaskl(i) = 0.0_kind_io8
           slmaskw(i) = 0.0_kind_io8
         ENDIF
       ENDIF

     ENDDO

   ELSE

! For running uncoupled (non-fractional grid).

     DO i=1,lensfc
       IF(nint(slmask(i)) == 1) THEN
         slmaskl(i) = 1.0_kind_io8
         slmaskw(i) = 1.0_kind_io8
       ELSE
         slmaskl(i) = 0.0_kind_io8
         slmaskw(i) = 0.0_kind_io8
       ENDIF
     ENDDO  

   ENDIF set_mask

   DO i=1,lensfc
     if(nint(slmask(i)) == 0) then
       min_ice(i) = 0.15_kind_io8
     else
       min_ice(i) = 0.0_kind_io8
     endif
   ENDDO

   num_threads = num_parthds()
   print*
   print*,"CALL SFCCYCLE TO UPDATE SURFACE FIELDS."
   CALL sfccycle(lugb,lensfc,lsoil,sig1t,deltsfc,          &
               iy,im,id,ih,fh,rla,rlo,                   &
               slmaskl,slmaskw, orog, orog_uf, use_ufo, do_nsst,   &
               sihfcs,sicfcs,sitfcs,sndfcs,slcfcs,       &
               vmnfcs,vmxfcs,slpfcs,absfcs,              &
               tsffcs,swefcs,zorfcs,albfcs,tg3fcs,       &
               cnpfcs,smcfcs,stcfcs,slifcs,aisfcs,       &
               vegfcs,vetfcs,sotfcs,alffcs,              &
               cvfcs,cvbfcs,cvtfcs,myrank,num_threads, nlunit,        &
               sz_nml, input_nml_file,                   &
               min_ice, &
               ialb,isot,ivegsrc,tile_num,i_index,j_index)

   DEALLOCATE(slmaskl, slmaskw)
 ENDIF

!--------------------------------------------------------------------------------
! IF RUNNING WITH NSST, READ IN GSI FILE WITH THE UPDATED INCREMENTS (ON THE
! GAUSSIAN GRID), INTERPOLATE INCREMENTS TO THE CUBED-SPHERE TILE, AND PERFORM
! REQUIRED ADJUSTMENTS AND QC.
!--------------------------------------------------------------------------------

 IF (do_nsst) THEN
   IF (nst_file == "NULL") THEN
     print*
     print*,"NO GSI FILE.  ADJUST IFD FOR FORMER ICE POINTS."
     DO i = 1, lensfc
       IF (sicfcs_fg(i) > 0.0 .AND. sicfcs(i) == 0) THEN
         nsst%IFD(i) = 3.0
       ENDIF
     ENDDO
     nsst%TFINC = 0.0
   ELSE
     print*
     print*,"ADJUST TREF FROM GSI INCREMENT"
!
!    Get tf climatology at the time
!
     call get_tf_clm(rla,rlo,jdim,idim,iy,im,id,ih,tf_clm,tf_trd)
     tf_clm_tile(:) = reshape(tf_clm, (/lensfc/) )
     tf_trd_tile(:) = reshape(tf_trd, (/lensfc/) )
!
!    Get salinity climatology at the time
!
     call get_sal_clm(rla,rlo,jdim,idim,iy,im,id,ih,sal_clm)
     sal_clm_tile(:) = reshape(sal_clm, (/lensfc/) )
!
!    read tf analysis increment generated by GSI
!
     CALL read_gsi_data(nst_file, 'NST')
!
!    update foundation & surface temperature for NSST
!
     CALL adjust_nsst(rla,rlo,slifcs,slifcs_fg,tsffcs,sitfcs,sicfcs,sicfcs_fg,&
                    stcfcs,nsst,lensfc,lsoil,idim,jdim,zsea1,zsea2, &
                    tf_clm_tile,tf_trd_tile,sal_clm_tile,landfrac,frac_grid)
   ENDIF
 ENDIF

!--------------------------------------------------------------------------------
! READ IN AND APPLY LAND INCREMENTS FROM THE GSI
!--------------------------------------------------------------------------------

 IF (do_lndinc) THEN

    ! SNOW INCREMENTS
    ! do snow first, as temperature updates will use snow analaysis
    IF (do_sno_inc) THEN
    ! updates are made to snow depth only over land (and not-land ice).
    ! SWE is then updated from the snow depth analysis, using the model
    ! forecast density

       !--------------------------------------------------------------------------------
       ! read increments in
       !--------------------------------------------------------------------------------

       ! Only coded for DA on native model grid (would always be the case for cycling DA)
       CALL read_data(lsoil,lensfc,.false.,.true.,sndfcs=snd_inc)

       !--------------------------------------------------------------------------------
       ! add increments to state vars
       !--------------------------------------------------------------------------------

       ! store background states
       snd_bck = sndfcs
       swe_bck = swefcs

       CALL add_increment_snow(snd_inc,landinc_mask,lensfc,sndfcs)

       !--------------------------------------------------------------------------------
       ! make any necessary adjustments to dependent variables
       !--------------------------------------------------------------------------------

       CALL apply_land_da_adjustments_snd(lsm, lensfc, landinc_mask, swe_bck, snd_bck, &
                        sndfcs, swefcs)

    ENDIF

    ! SOIL INCREMENTS
    IF (do_soi_inc) THEN

       !--------------------------------------------------------------------------------
       ! re-calculate soilsnow mask if snow has been updated.
       !--------------------------------------------------------------------------------

        landinc_mask_fg = landinc_mask

        IF (do_sfccycle .OR. do_sno_inc)  THEN
            CALL calculate_landinc_mask(smcfcs(:,1),swefcs, vetfcs, lensfc, &
                                                        veg_type_landice, landinc_mask )
        ENDIF

       !--------------------------------------------------------------------------------
       ! read increments in
       !--------------------------------------------------------------------------------

        WRITE(rankch, '(I3.3)') (myrank+1)

        lnd_soi_file = trim(lnd_soi_file) // "." //  rankch

        INQUIRE(file=trim(lnd_soi_file), exist=file_exists)
        IF (.not. file_exists) then
            print *, 'FATAL ERROR: land increment update requested, but file does not exist: ', &
                    trim(lnd_soi_file)
            call mpi_abort(mpi_comm_world, 10, ierr)
        ENDIF

        CALL read_gsi_data(lnd_soi_file, 'LND', lsoil=lsoil)

        !--------------------------------------------------------------------------------
        ! add increments to state vars
        !--------------------------------------------------------------------------------
        ! when applying increments, will often need to adjust other land states in response
        ! to the changes made. Need to store bacground, apply the increments, then make
        ! secondart adjustments. When updating more than one state, be careful about the
        ! order if increments and secondary adjustments.

        ! store background states
        stc_bck = stcfcs
        smc_bck = smcfcs
        slc_bck = slcfcs

        ! below updates [STC/SMC/STC]FCS to hold the analysis
        CALL add_increment_soil(rla,rlo,stcfcs,smcfcs,slcfcs,stc_updated, slc_updated, &
                landinc_mask,landinc_mask_fg,lensfc,lsoil,idim,jdim,lsm,myrank)

        !--------------------------------------------------------------------------------
        ! make any necessary adjustments to dependent variables
        !--------------------------------------------------------------------------------


        CALL apply_land_da_adjustments_soil(lsm, isot, ivegsrc,lensfc, lsoil, &
            sotfcs, landinc_mask_fg, stc_bck, stcfcs, smcfcs, slcfcs, stc_updated, &
            slc_updated,zsoil)

   ENDIF ! soil increments

!--------------------------------------------------------------------------------
! clean up
!--------------------------------------------------------------------------------

   ! to do - save and write out  STC_INC? (soil temperature increments)
   IF(ALLOCATED(landinc_mask_fg)) DEALLOCATE(landinc_mask_fg)
   IF(ALLOCATED(landinc_mask)) DEALLOCATE(landinc_mask)
   IF(ALLOCATED(stc_bck)) DEALLOCATE(stc_bck)
   IF(ALLOCATED(smc_bck)) DEALLOCATE(smc_bck)
   IF(ALLOCATED(slc_bck)) DEALLOCATE(slc_bck)
   IF(ALLOCATED(snd_bck)) DEALLOCATE(snd_bck)
   IF(ALLOCATED(swe_bck)) DEALLOCATE(swe_bck)
   IF(ALLOCATED(snd_inc)) DEALLOCATE(snd_inc)

 ENDIF
!--------------------------------------------------------------------------------
! WRITE OUT UPDATED SURFACE DATA ON THE CUBED-SPHERE TILE.
!--------------------------------------------------------------------------------

 IF (lsm==lsm_noahmp) THEN

   CALL write_data(lensfc,idim,jdim,lsoil,do_nsst,nsst,vegfcs=vegfcs, &
                   slcfcs=slcfcs,smcfcs=smcfcs,stcfcs=stcfcs,&
                   sicfcs=sicfcs,sihfcs=sihfcs)

 ELSEIF (lsm==lsm_noah) THEN

   CALL write_data(lensfc,idim,jdim,lsoil, &
                   do_nsst,nsst,slifcs=slifcs,tsffcs=tsffcs,vegfcs=vegfcs, &
                   swefcs=swefcs,tg3fcs=tg3fcs,zorfcs=zorfcs, &
                   albfcs=albfcs,alffcs=alffcs,cnpfcs=cnpfcs, &
                   f10m=f10m,t2m=t2m,q2m=q2m,vetfcs=vetfcs, &
                   sotfcs=sotfcs,ustar=ustar,fmm=fmm,fhh=fhh, &
                   sicfcs=sicfcs,sihfcs=sihfcs,sitfcs=sitfcs,tprcp=tprcp, &
                   srflag=srflag,swdfcs=sndfcs,vmnfcs=vmnfcs, &
                   vmxfcs=vmxfcs,slpfcs=slpfcs,absfcs=absfcs, &
                   slcfcs=slcfcs,smcfcs=smcfcs,stcfcs=stcfcs)

 ENDIF

 IF (do_nsst) THEN
   DEALLOCATE(nsst%C_0)
   DEALLOCATE(nsst%C_D)
   DEALLOCATE(nsst%D_CONV)
   DEALLOCATE(nsst%DT_COOL)
   DEALLOCATE(nsst%IFD)
   DEALLOCATE(nsst%QRAIN)
   DEALLOCATE(nsst%TREF)
   DEALLOCATE(nsst%TFINC)
   DEALLOCATE(nsst%W_0)
   DEALLOCATE(nsst%W_D)
   DEALLOCATE(nsst%XS)
   DEALLOCATE(nsst%XT)
   DEALLOCATE(nsst%XTTS)
   DEALLOCATE(nsst%XU)
   DEALLOCATE(nsst%XV)
   DEALLOCATE(nsst%XZ)
   DEALLOCATE(nsst%XZTS)
   DEALLOCATE(nsst%Z_C)
   DEALLOCATE(nsst%ZM)
   DEALLOCATE(slifcs_fg)
   DEALLOCATE(sicfcs_fg)
 ENDIF

 RETURN

 END SUBROUTINE sfcdrv
 
 SUBROUTINE adjust_nsst(RLA,RLO,SLMSK_TILE,SLMSK_FG_TILE,SKINT_TILE,&
                        SICET_TILE,sice_tile,sice_fg_tile,SOILT_TILE,NSST, &
                        LENSFC,LSOIL,IDIM,JDIM,ZSEA1,ZSEA2, &
                        tf_clm_tile,tf_trd_tile,sal_clm_tile,LANDFRAC, &
                        FRAC_GRID)

 USE utils
 USE gdswzd_mod
 USE read_write_data, ONLY : idim_gaus, jdim_gaus, &
                             slmsk_gaus, dtref_gaus, &
                             nsst_data

 USE mpi

 IMPLICIT NONE

 INTEGER, INTENT(IN)      :: LENSFC, LSOIL, IDIM, JDIM

 LOGICAL, INTENT(IN)      :: FRAC_GRID

 REAL, INTENT(IN)         :: SLMSK_TILE(LENSFC), SLMSK_FG_TILE(LENSFC), LANDFRAC(LENSFC)
 real, intent(in)         :: tf_clm_tile(lensfc),tf_trd_tile(lensfc),sal_clm_tile(lensfc)
 REAL, INTENT(IN)         :: ZSEA1, ZSEA2,sice_tile(lensfc),sice_fg_tile(lensfc)
 REAL, INTENT(IN)         :: RLA(LENSFC), RLO(LENSFC)
 REAL, INTENT(INOUT)      :: SKINT_TILE(LENSFC)
 REAL, INTENT(INOUT)      :: SICET_TILE(LENSFC),SOILT_TILE(LENSFC,LSOIL)

 TYPE(NSST_DATA)          :: NSST

 REAL, PARAMETER          :: TMAX=313.0,tzero=273.16

 INTEGER                  :: IOPT, NRET, KGDS_GAUS(200)
 INTEGER                  :: IGAUS, JGAUS, IJ, II, JJ, III, JJJ, KRAD
 INTEGER                  :: ISTART, IEND, JSTART, JEND
!INTEGER                  :: MASK_TILE, MASK_FG_TILE
 INTEGER,allocatable      :: MASK_TILE(:),MASK_FG_TILE(:)
 INTEGER                  :: ITILE, JTILE
 INTEGER                  :: MAX_SEARCH, J, IERR
 INTEGER                  :: IGAUSP1, JGAUSP1
 integer                  :: nintp,nsearched,nice,nland
 integer                  :: nfill,nfill_tice,nfill_clm
 integer                  :: nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c

 INTEGER, ALLOCATABLE     :: ID1(:,:), ID2(:,:), JDC(:,:)

 LOGICAL                  :: IS_ICE
 
 real                     :: tfreez
 REAL                     :: TREF_SAVE,WSUM,tf_ice,tf_thaw
 REAL                     :: FILL, DTZM, GAUS_RES_KM, DTREF
 REAL, ALLOCATABLE        :: XPTS(:), YPTS(:), LATS(:), LONS(:)
 REAL, ALLOCATABLE        :: DUM2D(:,:), LATS_RAD(:), LONS_RAD(:)
 REAL, ALLOCATABLE        :: AGRID(:,:,:), S2C(:,:,:)

 kgds_gaus     = 0
 kgds_gaus(1)  = 4          ! OCT 6 - TYPE OF GRID (GAUSSIAN)
 kgds_gaus(2)  = idim_gaus  ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
 kgds_gaus(3)  = jdim_gaus
 kgds_gaus(4)  = 90000      ! OCT 11-13 - LAT OF ORIGIN
 kgds_gaus(5)  = 0          ! OCT 14-16 - LON OF ORIGIN
 kgds_gaus(6)  = 128        ! OCT 17 - RESOLUTION FLAG
 kgds_gaus(7)  = -90000     ! OCT 18-20 - LAT OF EXTREME POINT
 kgds_gaus(8)  = nint(-360000./float(idim_gaus))  ! OCT 21-23 - LON OF EXTREME POINT
 kgds_gaus(9)  = nint((360.0 / float(idim_gaus))*1000.0)
                            ! OCT 24-25 - LONGITUDE DIRECTION INCR.
 kgds_gaus(10) = jdim_gaus/2     ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
 kgds_gaus(12) = 255        ! OCT 29 - RESERVED
 kgds_gaus(20) = 255        ! OCT 5  - NOT USED, SET TO 255

 print*
 print*,'ADJUST NSST USING GSI INCREMENTS ON GAUSSIAN GRID'

!----------------------------------------------------------------------
! CALL GDSWZD TO COMPUTE THE LAT/LON OF EACH GSI GAUSSIAN GRID POINT.
!----------------------------------------------------------------------

 iopt = 0
 fill = -9999.
 ALLOCATE(xpts(idim_gaus*jdim_gaus))
 ALLOCATE(ypts(idim_gaus*jdim_gaus))
 ALLOCATE(lats(idim_gaus*jdim_gaus))
 ALLOCATE(lons(idim_gaus*jdim_gaus))
 xpts = fill
 ypts = fill
 lats = fill
 lons = fill

 CALL gdswzd(kgds_gaus,iopt,(idim_gaus*jdim_gaus),fill,xpts,ypts,lons,lats,nret)

 IF (nret /= (idim_gaus*jdim_gaus)) THEN
   print*,'FATAL ERROR: PROBLEM IN GDSWZD. STOP.'
   CALL mpi_abort(mpi_comm_world, 12, ierr)
 ENDIF

 DEALLOCATE (xpts, ypts)

 ALLOCATE(dum2d(idim_gaus,jdim_gaus))
 dum2d = reshape(lats, (/idim_gaus,jdim_gaus/) )
 DEALLOCATE(lats)

 ALLOCATE(lats_rad(jdim_gaus))
 DO j = 1, jdim_gaus
   lats_rad(j) = dum2d(1,jdim_gaus-j+1) * 3.1415926 / 180.0
 ENDDO

 dum2d = reshape(lons, (/idim_gaus,jdim_gaus/) )
 DEALLOCATE(lons)
 ALLOCATE(lons_rad(idim_gaus))
 lons_rad = dum2d(:,1) * 3.1415926 / 180.0

 DEALLOCATE(dum2d)

 ALLOCATE(agrid(idim,jdim,2))
 agrid(:,:,1) = reshape(rlo, (/idim,jdim/) )
 agrid(:,:,2) = reshape(rla, (/idim,jdim/) )
 agrid        = agrid * 3.1415926 / 180.0

 ALLOCATE(id1(idim,jdim))
 ALLOCATE(id2(idim,jdim))
 ALLOCATE(jdc(idim,jdim))
 ALLOCATE(s2c(idim,jdim,4))

!----------------------------------------------------------------------
! COMPUTE BILINEAR WEIGHTS FOR EACH MODEL POINT FROM THE NEAREST
! FOUR GSI/GAUSSIAN POINTS.  DOES NOT ACCOUNT FOR MASK.  THAT
! HAPPENS LATER.
!----------------------------------------------------------------------

 CALL remap_coef( 1, idim, 1, jdim, idim_gaus, jdim_gaus, &
                  lons_rad, lats_rad, id1, id2, jdc, s2c, agrid )

 DEALLOCATE(lons_rad, lats_rad, agrid)

!----------------------------------------------------------------------
! THE MAXIMUM DISTANCE TO SEARCH IS 500 KM. HOW MANY GAUSSIAN
! GRID LENGTHS IS THAT?
!----------------------------------------------------------------------

 gaus_res_km = 360.0 / idim_gaus * 111.0
 max_search  = ceiling(500.0/gaus_res_km)

 print*
 print*,'MAXIMUM SEARCH IS ',max_search, ' GAUSSIAN POINTS.'
 print*

!
! Initialize variables for counts statitics to be zeros
!
 nintp = 0
 nset_thaw = 0
 nset_thaw_s = 0
 nset_thaw_i = 0
 nset_thaw_c = 0
 nsearched = 0
 nfill = 0
 nfill_tice = 0
 nfill_clm = 0
 nice = 0
 nland = 0
!----------------------------------------------------------------------
! TREF INCREMENT WILL BE OUTPUT.  INITIALIZE TO ZERO.
!----------------------------------------------------------------------

 nsst%TFINC = 0.0

 allocate(mask_tile(lensfc))
 allocate(mask_fg_tile(lensfc))
 
 IF(.NOT. frac_grid) THEN
   mask_tile    = nint(slmsk_tile)
   mask_fg_tile = nint(slmsk_fg_tile)
 ELSE
   mask_tile=0
   WHERE(sice_tile > 0.0) mask_tile=2
   WHERE(landfrac == 1.0) mask_tile=1
   mask_fg_tile=0
   WHERE(sice_fg_tile > 0.0) mask_fg_tile=2
   WHERE(landfrac == 1.0) mask_fg_tile=1
 ENDIF

 ij_loop : DO ij = 1, lensfc

!
!  when sea ice exists, get salinity dependent water temperature
!
   tf_ice = tfreez(sal_clm_tile(ij)) + tzero
!----------------------------------------------------------------------
! SKIP LAND POINTS.  NSST NOT APPLIED AT LAND.
!----------------------------------------------------------------------

   IF (mask_tile(ij) == 1) THEN
     nland = nland + 1
     cycle ij_loop  
   ENDIF

!
! these are ice points.  set tref to tf_ice and update tmpsfc.
!
   if (mask_tile(ij) == 2) then
     nsst%tref(ij)=tf_ice      ! water part tmp set
     skint_tile(ij)=(1.0-sice_tile(ij))*nsst%tref(ij)+sice_tile(ij)*sicet_tile(ij)
     nice = nice + 1
     cycle ij_loop
   endif

!
!  Get i,j index on array of (idim,jdim) from known ij
!
   jtile = (ij-1) / idim + 1
   itile = mod(ij,idim)
   IF (itile==0) itile = idim

!----------------------------------------------------------------------
! IF THE MODEL POINT WAS ICE COVERED, BUT IS NOW OPEN WATER, SET
! TREF TO searched adjascent open water onea, if failed the search, set to
! weighted average of tf_ice and tf_clm. For NSST vars, set xz TO '30' AND ALL OTHER FIELDS TO ZERO.
!----------------------------------------------------------------------

   IF (mask_fg_tile(ij) == 2 .AND. mask_tile(ij) == 0) THEN
!
!    set background for the thaw (just melted water) situation
!
     call tf_thaw_set(nsst%tref,mask_fg_tile,itile,jtile,tf_ice,tf_clm_tile(ij),tf_thaw,idim,jdim, &
                      nset_thaw_s,nset_thaw_i,nset_thaw_c)
     call nsst_water_reset(nsst,ij,tf_thaw)
     nset_thaw = nset_thaw + 1
   ENDIF

!----------------------------------------------------------------------
! THESE ARE POINTS THAT ARE OPEN WATER AND WERE OPEN WATER PRIOR
! TO ANY ICE UPDATE BY SFCCYCLE. UPDATE TREF AND SKIN TEMP.
! AT OPEN WATER POINTS, THE SEA ICE TEMPERATURE (SICET_TILE) AND
! SOIL COLUMN TEMPERATURE (SOILT_TILE) ARE SET TO THE SKIN TEMP.
! IT IS SIMPLY A FILLER VALUE.  THESE FIELDS ARE NOT USED AT
! OPEN WATER POINTS.
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! SEE IF ANY OF THE NEAREST GSI POINTS MASK AREA OPEN WATER.  
! IF SO, APPLY NSST INCREMENT USING BILINEAR INTERPOLATION.
!----------------------------------------------------------------------

   igaus   = id1(itile,jtile)
   jgaus   = jdc(itile,jtile)
   igausp1 = id2(itile,jtile)
   jgausp1 = jdc(itile,jtile)+1

   IF (slmsk_gaus(igaus,jgaus)     == 0 .OR. &
       slmsk_gaus(igausp1,jgaus)   == 0 .OR. &
       slmsk_gaus(igausp1,jgausp1) == 0 .OR. &
       slmsk_gaus(igaus,jgausp1)   == 0) THEN

     dtref = 0.0
     wsum  = 0.0

     IF (slmsk_gaus(igaus,jgaus) == 0) THEN
       dtref = dtref + (s2c(itile,jtile,1) * dtref_gaus(igaus,jgaus))
       wsum  = wsum + s2c(itile,jtile,1)
     ENDIF

     IF (slmsk_gaus(igausp1,jgaus) == 0) THEN
       dtref = dtref + (s2c(itile,jtile,2) * dtref_gaus(igausp1,jgaus))
       wsum  = wsum + s2c(itile,jtile,2)
     ENDIF

     IF (slmsk_gaus(igausp1,jgausp1) == 0) THEN
       dtref = dtref + (s2c(itile,jtile,3) * dtref_gaus(igausp1,jgausp1))
       wsum  = wsum + s2c(itile,jtile,3)
     ENDIF

     IF (slmsk_gaus(igaus,jgausp1) == 0) THEN
       dtref = dtref + (s2c(itile,jtile,4) * dtref_gaus(igaus,jgausp1))
       wsum  = wsum + s2c(itile,jtile,4)
     ENDIF

     nintp = nintp + 1
     dtref = dtref / wsum

     tref_save      = nsst%TREF(ij)
     nsst%TREF(ij)  = nsst%TREF(ij) + dtref
     nsst%TREF(ij)  = max(nsst%TREF(ij), tf_ice)
     nsst%TREF(ij)  = min(nsst%TREF(ij), tmax)
     nsst%TFINC(ij) = nsst%TREF(ij) - tref_save

     CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij),  &
                     nsst%Z_C(ij),zsea1,zsea2,dtzm)

     skint_tile(ij) = nsst%TREF(ij) + dtzm
     skint_tile(ij) = max(skint_tile(ij), tf_ice)
     skint_tile(ij) = min(skint_tile(ij), tmax)

     sicet_tile(ij)   = skint_tile(ij)
! Under fractional grids, soilt is used at points with at
! least some land.
     IF(.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)

!----------------------------------------------------------------------
! NO NEARBY GSI/GAUSSIAN OPEN WATER POINTS. PERFORM A SPIRAL SEARCH TO
! FIND NEAREST NON-LAND POINT ON GSI/GAUSSIAN GRID.
!----------------------------------------------------------------------

   ELSE  

     is_ice = .false.

     DO krad = 1, max_search

       istart = igaus - krad
       iend   = igaus + krad
       jstart = jgaus - krad
       jend   = jgaus + krad

       DO jj = jstart, jend
       DO ii = istart, iend

         IF((jj == jstart) .OR. (jj == jend) .OR.   &
            (ii == istart) .OR. (ii == iend))  THEN

           IF ((jj >= 1) .AND. (jj <= jdim_gaus)) THEN

             jjj = jj
             IF (ii <= 0) THEN
               iii = idim_gaus + ii
             ELSE IF (ii >= (idim_gaus+1)) THEN
               iii = ii - idim_gaus
             ELSE
               iii = ii
             END IF

!----------------------------------------------------------------------
! SEE IF NEARBY POINTS ARE SEA ICE.  IF THEY ARE, AND THE SEARCH FOR
! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
! FREEZING BELOW.
!----------------------------------------------------------------------

             IF (krad <= 2 .AND. slmsk_gaus(iii,jjj) == 2) is_ice = .true.

             IF (slmsk_gaus(iii,jjj) == 0) THEN

!              PRINT*,'MISMATCH AT TILE POINT  ',ITILE,JTILE
!              PRINT*,'UPDATE TREF USING GSI INCREMENT AT ',III,JJJ,DTREF_GAUS(III,JJJ)
               nsearched = nsearched + 1

               tref_save      = nsst%TREF(ij)
               nsst%TREF(ij ) = nsst%TREF(ij) + dtref_gaus(iii,jjj)
               nsst%TREF(ij)  = max(nsst%TREF(ij), tf_ice)
               nsst%TREF(ij)  = min(nsst%TREF(ij), tmax)
               nsst%TFINC(ij) = nsst%TREF(ij) - tref_save

               CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij),  &
                     nsst%Z_C(ij),zsea1,zsea2,dtzm)

               skint_tile(ij) = nsst%TREF(ij) + dtzm
               skint_tile(ij) = max(skint_tile(ij), tf_ice)
               skint_tile(ij) = min(skint_tile(ij), tmax)

               sicet_tile(ij)   = skint_tile(ij)
               IF(.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
               cycle ij_loop

             ENDIF ! GSI/Gaussian mask is open water

           ENDIF

         ENDIF

       ENDDO
       ENDDO

     ENDDO ! KRAD LOOP

!----------------------------------------------------------------------
! THE SEARCH FAILED.  IF THERE IS NEARBY ICE, SET TREF TO FREEZING.
! ELSE UPDATE TREF BASED ON THE ANNUAL SST CYCLE.
!----------------------------------------------------------------------

!    PRINT*,'WARNING !!!!!! SEARCH FAILED AT TILE POINT ',ITILE,JTILE

     nfill = nfill  + 1
     IF (is_ice) THEN
       nsst%TREF(ij) = tf_ice
!      PRINT*,"NEARBY ICE.  SET TREF TO FREEZING"
       nfill_tice = nfill_tice + 1
     ELSE
       tref_save      = nsst%TREF(ij)
       nsst%TREF(ij)  = nsst%TREF(ij) + tf_trd_tile(ij)
       nsst%TREF(ij)  = max(nsst%TREF(ij), tf_ice)
       nsst%TREF(ij)  = min(nsst%TREF(ij), tmax)
       nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
!      PRINT*,'UPDATE TREF FROM SST CLIMO ',DTREF
       nfill_clm = nfill_clm + 1
     ENDIF

     CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij),  &
                     nsst%Z_C(ij),zsea1,zsea2,dtzm)

     skint_tile(ij) = nsst%TREF(ij) + dtzm
     skint_tile(ij) = max(skint_tile(ij), tf_ice)
     skint_tile(ij) = min(skint_tile(ij), tmax)

     sicet_tile(ij)   = skint_tile(ij)
     IF (.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)

   ENDIF  ! NEARBY GAUSSIAN POINTS ARE OPEN WATER?

 ENDDO ij_loop

 write(*,'(a)') 'statistics of grids number processed for tile : '
 write(*,'(a,I8)') ' nintp = ',nintp
 write(*,'(a,4I8)') 'nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c =',nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c
 write(*,'(a,I8)') ' nsearched = ',nsearched
 write(*,'(a,3I6)') ' nfill,nfill_tice,nfill_clm = ',nfill,nfill_tice,nfill_clm
 write(*,'(a,I8)') ' nice = ',nice
 write(*,'(a,I8)') ' nland = ',nland

 DEALLOCATE(id1, id2, jdc, s2c, mask_tile, mask_fg_tile)

 END SUBROUTINE adjust_nsst

 SUBROUTINE climo_trend(LATITUDE, MON, DAY, DELTSFC, DTREF)
 IMPLICIT NONE

 INTEGER, INTENT(IN)    :: MON, DAY

 REAL, INTENT(IN)       :: LATITUDE, DELTSFC
 REAL, INTENT(OUT)      :: DTREF

 INTEGER                :: NUM_DAYS(12), MON2, MON1

 REAL, TARGET           :: SST_80_90(12)
 REAL, TARGET           :: SST_70_80(12)
 REAL, TARGET           :: SST_60_70(12)
 REAL, TARGET           :: SST_50_60(12)
 REAL, TARGET           :: SST_40_50(12)
 REAL, TARGET           :: SST_30_40(12)
 REAL, TARGET           :: SST_20_30(12)
 REAL, TARGET           :: SST_10_20(12)
 REAL, TARGET           :: SST_00_10(12)
 REAL, TARGET           :: SST_M10_00(12)
 REAL, TARGET           :: SST_M20_M10(12)
 REAL, TARGET           :: SST_M30_M20(12)
 REAL, TARGET           :: SST_M40_M30(12)
 REAL, TARGET           :: SST_M50_M40(12)
 REAL, TARGET           :: SST_M60_M50(12)
 REAL, TARGET           :: SST_M70_M60(12)
 REAL, TARGET           :: SST_M80_M70(12)
 REAL, TARGET           :: SST_M90_M80(12)

 REAL, POINTER          :: SST(:)

 DATA num_days  /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/

 DATA sst_80_90 /271.466, 271.458, 271.448, 271.445, 271.519, 271.636, & 
                 272.023, 272.066, 272.001, 271.698, 271.510, 271.472/

 DATA sst_70_80 /272.149, 272.103, 272.095, 272.126, 272.360, 272.988, &
                 274.061, 274.868, 274.415, 273.201, 272.468, 272.268/

 DATA sst_60_70 /274.240, 274.019, 273.988, 274.185, 275.104, 276.875, &
                 279.005, 280.172, 279.396, 277.586, 275.818, 274.803/

 DATA sst_50_60 /277.277, 276.935, 277.021, 277.531, 279.100, 281.357, &
                 283.735, 285.171, 284.399, 282.328, 279.918, 278.199/

 DATA sst_40_50 /281.321, 280.721, 280.850, 281.820, 283.958, 286.588, &
                 289.195, 290.873, 290.014, 287.652, 284.898, 282.735/

 DATA sst_30_40 /289.189, 288.519, 288.687, 289.648, 291.547, 293.904, &
                 296.110, 297.319, 296.816, 295.225, 292.908, 290.743/

 DATA sst_20_30 /294.807, 294.348, 294.710, 295.714, 297.224, 298.703, &
                 299.682, 300.127, 300.099, 299.455, 297.953, 296.177/

 DATA sst_10_20 /298.878, 298.720, 299.033, 299.707, 300.431, 300.709, &
                 300.814, 300.976, 301.174, 301.145, 300.587, 299.694/

 DATA sst_00_10 /300.415, 300.548, 300.939, 301.365, 301.505, 301.141, &
                 300.779, 300.660, 300.818, 300.994, 300.941, 300.675/

 DATA sst_m10_00 /300.226, 300.558, 300.914, 301.047, 300.645, 299.870, &
                  299.114, 298.751, 298.875, 299.294, 299.721, 299.989/

 DATA sst_m20_m10 /299.547, 299.985, 300.056, 299.676, 298.841, 297.788, &
                   296.893, 296.491, 296.687, 297.355, 298.220, 298.964/

 DATA sst_m30_m20 /297.524, 298.073, 297.897, 297.088, 295.846, 294.520, &
                   293.525, 293.087, 293.217, 293.951, 295.047, 296.363/

 DATA sst_m40_m30 /293.054, 293.765, 293.468, 292.447, 291.128, 289.781, &
                   288.773, 288.239, 288.203, 288.794, 289.947, 291.553/

 DATA sst_m50_m40 /285.052, 285.599, 285.426, 284.681, 283.761, 282.826, &
                   282.138, 281.730, 281.659, 281.965, 282.768, 283.961/

 DATA sst_m60_m50 /277.818, 278.174, 277.991, 277.455, 276.824, 276.229, &
                   275.817, 275.585, 275.560, 275.687, 276.142, 276.968/

 DATA sst_m70_m60 /273.436, 273.793, 273.451, 272.813, 272.349, 272.048, &
                   271.901, 271.838, 271.845, 271.889, 272.080, 272.607/

 DATA sst_m80_m70 /271.579, 271.578, 271.471, 271.407, 271.392, 271.391, &
                   271.390, 271.391, 271.394, 271.401, 271.422, 271.486/

 DATA sst_m90_m80 /271.350, 271.350, 271.350, 271.350, 271.350, 271.350, &
                   271.350, 271.350, 271.350, 271.350, 271.350, 271.350/

 NULLIFY(sst)
 IF (latitude > 80.0) THEN
   sst => sst_80_90
 ELSEIF (latitude > 70.0) THEN
   sst => sst_70_80
 ELSEIF (latitude > 60.0) THEN
   sst => sst_60_70
 ELSEIF (latitude > 50.0) THEN
   sst => sst_50_60
 ELSEIF (latitude > 40.0) THEN
   sst => sst_40_50
 ELSEIF (latitude > 30.0) THEN
   sst => sst_30_40
 ELSEIF (latitude > 20.0) THEN
   sst => sst_20_30
 ELSEIF (latitude > 10.0) THEN
   sst => sst_10_20
 ELSEIF (latitude > 0.0) THEN
   sst => sst_00_10
 ELSEIF (latitude > -10.0) THEN
   sst => sst_m10_00
 ELSEIF (latitude > -20.0) THEN
   sst => sst_m20_m10
 ELSEIF (latitude > -30.0) THEN
   sst => sst_m30_m20
 ELSEIF (latitude > -40.0) THEN
   sst => sst_m40_m30
 ELSEIF (latitude > -50.0) THEN
   sst => sst_m50_m40
 ELSEIF (latitude > -60.0) THEN
   sst => sst_m60_m50
 ELSEIF (latitude > -70.0) THEN
   sst => sst_m70_m60
 ELSEIF (latitude > -80.0) THEN
   sst => sst_m80_m70
 ELSE
   sst => sst_m90_m80
 END IF

 IF (day >= 15) THEN
   mon2 = mon+1
   IF(mon2 == 13) mon2 = 1
   mon1 = mon
   dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
 ELSE
   mon1 = mon - 1
   IF (mon1 == 0) mon1=12
   mon2 = mon
   dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
 ENDIF

 dtref = dtref * (deltsfc / 24.0)

 END SUBROUTINE climo_trend

 SUBROUTINE dtzm_point(XT,XZ,DT_COOL,ZC,Z1,Z2,DTZM)
  implicit none

  real, intent(in)  :: xt,xz,dt_cool,zc,z1,z2
  real, intent(out) :: dtzm

  real, parameter   :: zero = 0.0
  real, parameter   :: one  = 1.0
  real, parameter   :: half = 0.5
  real              :: dt_warm,dtw,dtc
!
! get the mean warming in the range of z=z1 to z=z2
!
  dtw = zero
  if ( xt > zero ) then
    dt_warm = (xt+xt)/xz      ! Tw(0)
    if ( z1 < z2) then
      if ( z2 < xz ) then
        dtw = dt_warm*(one-(z1+z2)/(xz+xz))
      elseif ( z1 < xz .and. z2 >= xz ) then
        dtw = half*(one-z1/xz)*dt_warm*(xz-z1)/(z2-z1)
      endif
    elseif ( z1 == z2 ) then
      if ( z1 < xz ) then
        dtw = dt_warm*(one-z1/xz)
      endif
    endif
  endif
!
! get the mean cooling in the range of z=z1 to z=z2
!
  dtc = zero
  if ( zc > zero ) then
    if ( z1 < z2) then
      if ( z2 < zc ) then
        dtc = dt_cool*(one-(z1+z2)/(zc+zc))
      elseif ( z1 < zc .and. z2 >= zc ) then
        dtc = half*(one-z1/zc)*dt_cool*(zc-z1)/(z2-z1)
      endif
    elseif ( z1 == z2 ) then
      if ( z1 < zc ) then
        dtc = dt_cool*(one-z1/zc)
      endif
    endif
  endif

!
! get the mean T departure from Tf in the range of z=z1 to z=z2
!
  dtzm = dtw - dtc

 END SUBROUTINE dtzm_point

 subroutine tf_thaw_set(tf_ij,mask_ij,itile,jtile,tice,tclm,tf_thaw,nx,ny, &
                        nset_thaw_s,nset_thaw_i,nset_thaw_c)

 real,    dimension(nx*ny), intent(in)    :: tf_ij
 integer, dimension(nx*ny), intent(in)    :: mask_ij
 real,                      intent(in)    :: tice,tclm
 integer,                   intent(in)    :: itile,jtile,nx,ny
 real,                      intent(out)   :: tf_thaw
 integer,                   intent(inout) :: nset_thaw_s,nset_thaw_i,nset_thaw_c
! Local
 real, parameter :: bmiss = -999.0
 real,    dimension(nx,ny) :: tf
 integer, dimension(nx,ny) :: mask
 integer :: krad,max_search,istart,iend,jstart,jend
 integer :: ii,jj,iii,jjj
 logical :: is_ice

 max_search = 2

 mask(:,:) = reshape(mask_ij,(/nx,ny/) )
 tf(:,:)   = reshape(tf_ij,(/nx,ny/) )

 tf_thaw = bmiss

 do krad = 1, max_search

    istart = itile - krad
    iend   = itile + krad
    jstart = jtile - krad
    jend   = jtile + krad

    do jj = jstart, jend
       do ii = istart, iend


          if ((jj == jstart) .or. (jj == jend) .or.   &
              (ii == istart) .or. (ii == iend))  then

             if ((jj >= 1) .and. (jj <= ny)) then
                jjj = jj
                if (ii <= 0) then
                   iii = nx + ii
                else if (ii >= (nx+1)) then
                   iii = ii - nx
                else
                   iii = ii
                endif

!----------------------------------------------------------------------
! SEE IF NEARBY POINTS ARE SEA ICE.  IF THEY ARE, AND THE SEARCH FOR
! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
! FREEZING BELOW.
!----------------------------------------------------------------------
                if (krad <= 2 .and. mask(iii,jjj) == 2) is_ice = .true.

                if (mask(iii,jjj) == 0) then
                   tf_thaw = tf(iii,jjj)
                   nset_thaw_s = nset_thaw_s + 1
                   write(*,'(a,I4,2F9.3)') 'nset_thaw_s,tf(iii,jjj),tclm : ',nset_thaw_s,tf(iii,jjj),tclm
                   go to 100
                endif ! tile mask is open water

             endif
          endif
       enddo
    enddo
 enddo  ! krad loop

 100 continue

 if ( tf_thaw == bmiss ) then
    if (is_ice) then
       tf_thaw = tice
       nset_thaw_i = nset_thaw_i + 1
       write(*,'(a,I4,F9.3)') 'nset_thaw_i,tf_ice : ',nset_thaw_i,tice
    else
       tf_thaw = 0.8*tice+0.2*tclm
       nset_thaw_c = nset_thaw_c + 1
       write(*,'(a,I4,2F9.3)') 'nset_thaw_c,tf_ice,tclm : ',nset_thaw_c,tice,tclm
    endif
 endif

 end subroutine tf_thaw_set

 subroutine nsst_water_reset(nsst,ij,tf_thaw)
 use read_write_data, only : nsst_data
 implicit none

 integer, intent(in)  :: ij

 real, intent(in)     :: tf_thaw

 type(nsst_data), intent(inout)      :: nsst

 nsst%c_0(ij)     = 0.0
 nsst%c_d(ij)     = 0.0
 nsst%d_conv(ij)  = 0.0
 nsst%dt_cool(ij) = 0.0
 nsst%ifd(ij)     = 0.0
 nsst%qrain(ij)   = 0.0
 nsst%tref(ij)    = tf_thaw
 nsst%w_0(ij)     = 0.0
 nsst%w_d(ij)     = 0.0
 nsst%xs(ij)      = 0.0
 nsst%xt(ij)      = 0.0
 nsst%xtts(ij)    = 0.0
 nsst%xu(ij)      = 0.0
 nsst%xv(ij)      = 0.0
 nsst%xz(ij)      = 30.0
 nsst%xzts(ij)    = 0.0
 nsst%z_c(ij)     = 0.0
 nsst%zm(ij)      = 0.0

 end subroutine nsst_water_reset

subroutine get_tf_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,tf_clm,tf_trd)
 use read_write_data, only : get_tf_clm_dim

 implicit none

 real,    dimension(nx*ny), intent(in)  :: xlats_ij 
 real,    dimension(nx*ny), intent(in)  :: xlons_ij 
 real,    dimension(nx,ny), intent(out) :: tf_clm   
 real,    dimension(nx,ny), intent(out) :: tf_trd   
 integer, intent(in) :: iy,im,id,ih,nx,ny
! local declare
 real,    allocatable, dimension(:,:)   :: tf_clm0    ! sst climatology at the valid time (nxc,nyc)
 real,    allocatable, dimension(:,:)   :: tf_trd0    ! 6-hourly sst climatology tendency valid at atime (nxc,nyc)
 real,    allocatable, dimension(:)     :: cxlats     ! latitudes of sst climatology
 real,    allocatable, dimension(:)     :: cxlons     ! longitudes of sst climatology

 real,    dimension(nx*ny)  :: tf_clm_ij  ! sst climatology at target grids (nx*ny)
 real,    dimension(nx*ny)  :: tf_trd_ij  ! 6-hourly sst climatology tendency 
 real :: wei1,wei2
 integer :: nxc,nyc,mon1,mon2
 character (len=6), parameter :: fin_tf_clm='sstclm' ! sst climatology file name
!
! get which two months used and their weights from atime
!
 call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
!
! get the dimensions of the sst climatology & allocate the related arrays
!
 call get_tf_clm_dim(fin_tf_clm,nyc,nxc)
 allocate( tf_clm0(nxc,nyc),tf_trd0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
!
! get tf_clm at the analysis time from monthly climatology & cxlats, cxlons
!
 call get_tf_clm_ta(tf_clm0,tf_trd0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
!
! get tf_clm (nx by ny lat/lon) valid at atime
!
 if ( nx == nxc .and. ny == nyc ) then
    tf_clm(:,:) = tf_clm0(:,:)
    tf_trd(:,:) = tf_trd0(:,:)
!   write(*,'(a,2f9.3)') 'same dimensions, tf_clm, min : ',minval(tf_clm),maxval(tf_clm)
 else
!   write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
    call intp_tile(tf_clm0,  cxlats,  cxlons,  nyc, nxc, &
                   tf_clm_ij,xlats_ij,xlons_ij,ny,  nx)
    call intp_tile(tf_trd0,  cxlats,  cxlons,  nyc, nxc, &
                   tf_trd_ij,xlats_ij,xlons_ij,ny,  nx)
!   write(*,'(a,2f9.3)') 'tf_clm0, min, max                        : ',minval(tf_clm0),maxval(tf_clm0)

    tf_clm(:,:) = reshape(tf_clm_ij, (/nx,ny/) )
    tf_trd(:,:) = reshape(tf_trd_ij, (/nx,ny/) )
 endif

end subroutine get_tf_clm

subroutine get_tf_clm_ta(tf_clm_ta,tf_clm_trend,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)
 use read_write_data, only : read_tf_clim_grb
 implicit none

! input
 integer, intent(in) :: nlat,nlon,mon1,mon2
 real,    intent(in) :: wei1,wei2
! output
 real, dimension(nlon,nlat), intent(out)   :: tf_clm_ta,tf_clm_trend
 real, dimension(nlat),      intent(out)   :: xlats
 real, dimension(nlon),      intent(out)   :: xlons

!input/output data file names
 character (len=6),  parameter :: fin_tf_clm='sstclm'

! local declare
 real, dimension(nlon,nlat) :: tf_clm1,tf_clm2

!
! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
!
  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm1,xlats,xlons,nlat,nlon,mon1)
  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm2,xlats,xlons,nlat,nlon,mon2)
!
!  tf_clim at the analysis time
!
   tf_clm_ta(:,:) = wei1*tf_clm1(:,:)+wei2*tf_clm2(:,:)
!
!  tf_clim trend at the analysis time (month to month)
!
   tf_clm_trend(:,:) = (tf_clm2(:,:)-tf_clm1(:,:))/120.0

   write(*,'(a,2f9.3)') 'tf_clm_ta, min, max : ',minval(tf_clm_ta),maxval(tf_clm_ta)
   write(*,'(a,2f9.3)') 'tf_clm_trend, min, max : ',minval(tf_clm_trend),maxval(tf_clm_trend)
 end subroutine get_tf_clm_ta

subroutine get_sal_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,sal_clm)
 use read_write_data, only : get_dim_nc
 implicit none

 real,    dimension(nx*ny), intent(in)  :: xlats_ij   ! 
 real,    dimension(nx*ny), intent(in)  :: xlons_ij   ! 
 real,    dimension(nx,ny), intent(out) :: sal_clm    ! 
 integer, intent(in) :: iy,im,id,ih,nx,ny
! local declare
 real,    allocatable, dimension(:,:)   :: sal_clm0   ! salinity climatology at the valid time
 real,    allocatable, dimension(:)     :: cxlats     ! latitudes of sst climatology
 real,    allocatable, dimension(:)     :: cxlons     ! longitudes of sst climatology

 real,    dimension(nx*ny)  :: sal_clm_ij  ! salinity climatology at target grids (nx*ny)
 real :: wei1,wei2
 integer :: nxc,nyc,mon1,mon2
 character (len=6), parameter :: fin_sal_clm='salclm' ! salinity climatology file name
!
! get which two months used and their weights from atime
!
 call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
!
! get the dimensions of the sst climatology & allocate the related arrays
!
 call get_dim_nc(fin_sal_clm,nyc,nxc)
 allocate( sal_clm0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
!
! get sal_clm at the analysis time from monthly climatology & cxlats, cxlons
!
 call get_sal_clm_ta(sal_clm0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
!
! get sal_clm (nx by ny lat/lon) valid at atime
!
 if ( nx == nxc .and. ny == nyc ) then
    sal_clm(:,:) = sal_clm0(:,:)
!   write(*,'(a,2f9.3)') 'same dimensions, sal_clm, min : ',minval(sal_clm),maxval(sal_clm)
 else
!   write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
    call intp_tile(sal_clm0,  cxlats,  cxlons,  nyc, nxc, &
                   sal_clm_ij,xlats_ij,xlons_ij,ny,  nx)
!   write(*,'(a,2f9.3)') 'sal_clm0, min, max                        : ',minval(sal_clm0),maxval(sal_clm0)
!   write(*,'(a,2f9.3)') 'done with intp_tile for sal_clm, min, max : ',minval(sal_clm_ij),maxval(sal_clm_ij)

    sal_clm(:,:) = reshape(sal_clm_ij, (/nx,ny/) )
 endif

end subroutine get_sal_clm

subroutine get_sal_clm_ta(sal_clm_ta,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)

 use read_write_data, only : read_salclm_gfs_nc
 implicit none

! input
 integer, intent(in) :: nlat,nlon,mon1,mon2
 real,    intent(in) :: wei1,wei2
! output
 real, dimension(nlon,nlat), intent(out)   :: sal_clm_ta
 real, dimension(nlat),      intent(out)   :: xlats
 real, dimension(nlon),      intent(out)   :: xlons

!input/output data file names
 character (len=6),  parameter :: fin_sal_clm='salclm'

! local declare
 real, dimension(nlon,nlat) :: sal_clm1,sal_clm2

!
! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
!
  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm1,xlats,xlons,nlat,nlon,mon1)
  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm2,xlats,xlons,nlat,nlon,mon2)
!
!  sal_clim at the analysis time
!
   sal_clm_ta(:,:) = wei1*sal_clm1(:,:)+wei2*sal_clm2(:,:)
   write(*,'(a,2f9.3)') 'sal_clm_ta, min, max : ',minval(sal_clm_ta),maxval(sal_clm_ta)
 end subroutine get_sal_clm_ta

subroutine intp_tile(tf_lalo,dlats_lalo,dlons_lalo,jdim_lalo,idim_lalo, &
                     tf_tile,xlats_tile,xlons_tile,jdim_tile,idim_tile)

 use utils

 implicit none

! input/output
 real,    dimension(idim_lalo,jdim_lalo), intent(in)  :: tf_lalo
 real,    dimension(jdim_lalo),           intent(in)  :: dlats_lalo
 real,    dimension(idim_lalo),           intent(in)  :: dlons_lalo
 real,    dimension(jdim_tile*idim_tile), intent(in)  :: xlats_tile
 real,    dimension(jdim_tile*idim_tile), intent(in)  :: xlons_tile
 integer,                                 intent(in)  :: jdim_lalo,idim_lalo,jdim_tile,idim_tile
 real,    dimension(jdim_tile*idim_tile), intent(out) :: tf_tile

! local
 real, parameter :: deg2rad=3.1415926/180.0
 real,    dimension(jdim_lalo) :: xlats_lalo
 real,    dimension(idim_lalo) :: xlons_lalo
 real    :: wsum
 integer :: itile,jtile
 integer :: ij
 integer :: ilalo,jlalo,ilalop1,jlalop1

 integer, allocatable, dimension(:,:)   :: id1,id2,jdc
 real,    allocatable, dimension(:,:,:) :: agrid,s2c

 print*
 print*,'interpolate from lat/lon grids to any one grid with known lat/lon'

 xlats_lalo = dlats_lalo*deg2rad
 xlons_lalo = dlons_lalo*deg2rad

 allocate(agrid(idim_tile,jdim_tile,2))
 agrid(:,:,1) = reshape(xlons_tile, (/idim_tile,jdim_tile/) )
 agrid(:,:,2) = reshape(xlats_tile, (/idim_tile,jdim_tile/) )
 agrid        = agrid*deg2rad

 allocate(id1(idim_tile,jdim_tile))
 allocate(id2(idim_tile,jdim_tile))
 allocate(jdc(idim_tile,jdim_tile))
 allocate(s2c(idim_tile,jdim_tile,4))

!----------------------------------------------------------------------
! compute bilinear weights for each model point from the nearest
! four lalo points. does not account for mask.  that
! happens later.
!----------------------------------------------------------------------

 call remap_coef( 1, idim_tile, 1, jdim_tile, idim_lalo, jdim_lalo, &
                  xlons_lalo, xlats_lalo, id1, id2, jdc, s2c, agrid )

 do ij = 1, jdim_tile*idim_tile

    jtile = (ij-1)/idim_tile + 1
    itile = mod(ij,idim_tile)
    if (itile==0) itile = idim_tile

    ilalo   = id1(itile,jtile)
    ilalop1 = id2(itile,jtile)
    jlalo   = jdc(itile,jtile)
    jlalop1 = jdc(itile,jtile) + 1

    wsum = s2c(itile,jtile,1) + s2c(itile,jtile,2) + &
           s2c(itile,jtile,3) + s2c(itile,jtile,4)

    tf_tile(ij)  = ( s2c(itile,jtile,1)*tf_lalo(ilalo,jlalo)     + &
                     s2c(itile,jtile,2)*tf_lalo(ilalop1,jlalo)   + &
                     s2c(itile,jtile,3)*tf_lalo(ilalop1,jlalop1) + &
                     s2c(itile,jtile,4)*tf_lalo(ilalo,jlalop1) )/wsum
 enddo

 deallocate(id1, id2, jdc, s2c)

end subroutine intp_tile

subroutine get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
 implicit none

! input
 integer, intent(in) :: iy,im,id,ih
! output
 integer, intent(out) :: mon1,mon2
 real,    intent(out) :: wei1,wei2

! local declare
 real :: rjday
 integer :: mon,monend,monm,monp,jdow,jdoy,jday
 integer :: jda(8)
!
!dayhf : julian day of the middle of each month
!
 real, dimension(13) ::  dayhf
 data dayhf/15.5,45.0,74.5,105.0,135.5,166.0,196.5,227.5,258.0,288.5,319.0,349.5,380.5/

! 15, 44, 73.5, 104, 134.5, 165, 195.5, 226.5, 257, 287.5, 318.5, 349  ! from
! woa05

 jda=0
 jda(1)=iy
 jda(2)=im
 jda(3)=id
 jda(5)=ih
 jdow = 0
 jdoy = 0
 jday = 0
 call w3doxdat(jda,jdow,jdoy,jday)
 rjday=jdoy+jda(5)/24.
 if(rjday.lt.dayhf(1)) rjday=rjday+365.

 monend = 12
 do mon = 1, monend
    monm = mon
    monp = mon + 1
    if( rjday >= dayhf(monm) .and. rjday < dayhf(monp) ) then
       mon1 = monm
       mon2 = monp
       go to 10
    endif
 enddo

 print *,'FATAL ERROR in get_tim_wei, wrong rjday',rjday
 call abort
 10     continue

 wei1 = (dayhf(mon2)-rjday)/(dayhf(mon2)-dayhf(mon1))
 wei2 = (rjday-dayhf(mon1))/(dayhf(mon2)-dayhf(mon1))

 if( mon2 == 13 ) mon2=1

 write(*,'(a,2i4,3f9.3)') 'mon1,mon2,rjday,wei1,wei2=',mon1,mon2,rjday,wei1,wei2

 end subroutine get_tim_wei

 real function tfreez(salinity)

 implicit none

 REAL salinity,sal
 REAL a1, a2, a3
 parameter(a1 = -0.0575)
 parameter(a2 =  1.710523e-3)
 parameter(a3 = -2.154996e-4)

 IF (salinity .LT. 0.) THEN
   sal = 0.
 ELSE
   sal = salinity
 ENDIF
 tfreez = sal*(a1+a2*sqrt(sal)+a3*sal)

 return
 end