ecearth3/sources/nemo-3.6/NEMO/OOO_SRC/ooo_read.F90

MODULE ooo_read
   !!==================================================================
   !!                    *** MODULE ooo_read ***
   !! Read routines : I/O for offline obs_oper
   !!==================================================================

   USE mppini
   USE lib_mpp
   USE in_out_manager
   USE par_kind, ONLY: lc
   USE netcdf
   USE oce,     ONLY: tsn, sshn
   USE dom_oce, ONLY: nlci, nlcj, nimpp, njmpp, tmask
   USE par_oce, ONLY: jpi, jpj, jpk
   USE obs_fbm, ONLY: fbimdi, fbrmdi, fbsp, fbdp
   USE ooo_data

   IMPLICIT NONE
   PRIVATE

   PUBLIC ooo_rea_dri

   !! $Id: ooo_read.F90 2355 2015-05-20 07:11:50Z ufla $
CONTAINS
   SUBROUTINE ooo_rea_dri(kfile)
      !!------------------------------------------------------------------------
      !!             *** ooo_rea_dri ***
      !!
      !! Purpose : To choose appropriate read method
      !! Method  : 
      !!
      !! Author  : A. Ryan Oct 2013
      !!
      !!------------------------------------------------------------------------
      INTEGER, INTENT(IN) :: &
              & kfile         !: File number
      CHARACTER(len=lc) :: &
              & cdfilename, & !: File name
              & cmatchname    !: Match name
      INTEGER :: &
              & kindex       !: File index to read
 
      !! Filename, index and match-up kind
      cdfilename = TRIM(ooo_files(kfile))
      cmatchname = TRIM(cl4_vars(kfile))
      kindex = nn_ooo_idx(kfile)

      !! Update model fields
      !! Class 4 variables: forecast, persistence,
      !!                    nrt_analysis, best_estimate
      !! Feedback variables: empty string
      IF ( (TRIM(cmatchname) == 'forecast') .OR. &
         & (TRIM(cmatchname) == 'persistence') .OR. &
         & (TRIM(cmatchname) == 'nrt_analysis') .OR. &
         & (TRIM(cmatchname) == 'best_estimate').OR. &
         & (TRIM(cmatchname) == '') ) THEN
         CALL ooo_read_file(TRIM(cdfilename), kindex)
         CALL ooo_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
      ELSE IF (TRIM(cmatchname) == 'climatology') THEN
         WRITE(numout,*) 'Interpolating climatologies'
      ELSE IF (TRIM(cmatchname) == 'altimeter') THEN
         CALL ooo_read_altbias(TRIM(cdfilename))
         CALL ooo_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
      END IF

   END SUBROUTINE ooo_rea_dri

   SUBROUTINE ooo_read_altbias(filename)
      !!------------------------------------------------------------------------
      !!                      *** ooo_read_altbias ***
      !!
      !! Purpose : To read altimeter bias and set tn,sn to missing values
      !! Method  : Use subdomain indices to create start and count matrices
      !!           for netcdf read.
      !!
      !! Author  : A. Ryan Sept 2012
      !!------------------------------------------------------------------------
      CHARACTER(len=*), INTENT(IN) :: filename
      INTEGER                      :: ncid, &
                                    & varid,&
                                    & istat,&
                                    & ntimes,&
                                    & tdim, &
                                    & xdim, &
                                    & ydim, &
                                    & zdim
      INTEGER                      :: ii, ij, ik
      INTEGER, DIMENSION(3)        :: start_s, &
                                    & count_s
      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
      REAL(fbdp)                     :: fill_val

      IF (TRIM(filename) == 'nofile') THEN
         tsn(:,:,:,:) = fbrmdi
         sshn(:,:) = fbrmdi 
      ELSE
         ! Open Netcdf file to find dimension id
         istat = nf90_open(trim(filename),nf90_nowrite,ncid)
         istat = nf90_inq_dimid(ncid,'x',xdim)
         istat = nf90_inq_dimid(ncid,'y',ydim)
         istat = nf90_inq_dimid(ncid,'deptht',zdim)
         istat = nf90_inq_dimid(ncid,'time',tdim)
         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
         ! Allocate temporary temperature array
         ALLOCATE(temp_sshn(nlci,nlcj))
         ! Create start and count arrays
         start_s = (/ nimpp, njmpp, 1 /)
         count_s = (/ nlci,  nlcj,  1 /)
          
         ! Altimeter bias
         istat = nf90_inq_varid(ncid,'altbias',varid)
         istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
         istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s)
         WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi
   
         ! Initialise tsn, sshn to fbrmdi
         tsn(:,:,:,:) = fbrmdi
         sshn(:,:) = fbrmdi 

         ! Fill sshn with altimeter bias 
         sshn(1:nlci,1:nlcj) = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)

         ! Remove halo from tmask, sshn to prevent double obs counting
         IF (jpi > nlci) THEN
             tmask(nlci+1:,:,:) = 0
             sshn(nlci+1:,:) = 0
         END IF
         IF (jpj > nlcj) THEN
             tmask(:,nlcj+1:,:) = 0
             sshn(:,nlcj+1:) = 0
         END IF

         ! Deallocate arrays
         DEALLOCATE(temp_sshn)

         ! Close netcdf file
         istat = nf90_close(ncid)
      END IF
   
   END SUBROUTINE ooo_read_altbias

   SUBROUTINE ooo_read_file(filename, ifcst)
      !!------------------------------------------------------------------------
      !!             *** ooo_read_file ***
      !!
      !! Purpose : To fill tn and sn with dailymean field from netcdf files
      !! Method  : Use subdomain indices to create start and count matrices
      !!           for netcdf read.
      !!
      !! Author  : A. Ryan Oct 2010
      !!------------------------------------------------------------------------

      INTEGER,          INTENT(IN) :: ifcst
      CHARACTER(len=*), INTENT(IN) :: filename
      INTEGER                      :: ncid, &
                                    & varid,&
                                    & istat,&
                                    & ntimes,&
                                    & tdim, &
                                    & xdim, &
                                    & ydim, &
                                    & zdim
      INTEGER                      :: ii, ij, ik
      INTEGER, DIMENSION(4)        :: start_n, &
                                    & count_n
      INTEGER, DIMENSION(3)        :: start_s, &
                                    & count_s
      REAL(fbdp), DIMENSION(:,:,:),ALLOCATABLE :: temp_tn, &
                                              & temp_sn
      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
      REAL(fbdp)                     :: fill_val

      ! DEBUG
      INTEGER :: istage

      IF (TRIM(filename) == 'nofile') THEN
         tsn(:,:,:,:) = fbrmdi
         sshn(:,:) = fbrmdi 
      ELSE
         WRITE(numout,*) "Opening :", TRIM(filename)
         ! Open Netcdf file to find dimension id
         istat = nf90_open(path=TRIM(filename), mode=nf90_nowrite, ncid=ncid)
         IF ( istat /= nf90_noerr ) THEN
             WRITE(numout,*) "WARNING: Could not open ", trim(filename)
             WRITE(numout,*) "ERROR: ", nf90_strerror(istat)
         ENDIF
         istat = nf90_inq_dimid(ncid,'x',xdim)
         istat = nf90_inq_dimid(ncid,'y',ydim)
         istat = nf90_inq_dimid(ncid,'deptht',zdim)
         istat = nf90_inq_dimid(ncid,'time_counter',tdim)
         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
         IF (ifcst .LE. ntimes) THEN   
            ! Allocate temporary temperature array
            ALLOCATE(temp_tn(nlci,nlcj,jpk))
            ALLOCATE(temp_sn(nlci,nlcj,jpk))
            ALLOCATE(temp_sshn(nlci,nlcj))
      
            ! Set temp_tn, temp_sn to 0.
            temp_tn(:,:,:) = fbrmdi
            temp_sn(:,:,:) = fbrmdi
            temp_sshn(:,:) = fbrmdi
      
            ! Create start and count arrays
            start_n = (/ nimpp, njmpp, 1,   ifcst /)
            count_n = (/ nlci,  nlcj,  jpk, 1     /)
            start_s = (/ nimpp, njmpp, ifcst /)
            count_s = (/ nlci,  nlcj,  1     /)
             
            ! Read information into temporary arrays
            ! retrieve varid and read in temperature
            istat = nf90_inq_varid(ncid,'votemper',varid)
            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
            istat = nf90_get_var(ncid, varid, temp_tn, start_n, count_n)
            WHERE(temp_tn(:,:,:) == fill_val) temp_tn(:,:,:) = fbrmdi
      
            ! retrieve varid and read in salinity
            istat = nf90_inq_varid(ncid,'vosaline',varid)
            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
            istat = nf90_get_var(ncid, varid, temp_sn, start_n, count_n)
            WHERE(temp_sn(:,:,:) == fill_val) temp_sn(:,:,:) = fbrmdi
      
            ! retrieve varid and read in SSH
            istat = nf90_inq_varid(ncid,'sossheig',varid)
            IF (istat /= nf90_noerr) THEN
               ! Altimeter bias
               istat = nf90_inq_varid(ncid,'altbias',varid)
            END IF
      
            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
            istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s)
            WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi
   
            ! Initialise tsn, sshn to fbrmdi
            tsn(:,:,:,:) = fbrmdi
            sshn(:,:) = fbrmdi 

            ! Mask out missing data index
            tsn(1:nlci,1:nlcj,1:jpk,1) = temp_tn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
            tsn(1:nlci,1:nlcj,1:jpk,2) = temp_sn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
            sshn(1:nlci,1:nlcj)        = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)

            ! Remove halo from tmask, tsn, sshn to prevent double obs counting
            IF (jpi > nlci) THEN
                tmask(nlci+1:,:,:) = 0
                tsn(nlci+1:,:,:,1) = 0
                tsn(nlci+1:,:,:,2) = 0
                sshn(nlci+1:,:) = 0
            END IF
            IF (jpj > nlcj) THEN
                tmask(:,nlcj+1:,:) = 0
                tsn(:,nlcj+1:,:,1) = 0
                tsn(:,nlcj+1:,:,2) = 0
                sshn(:,nlcj+1:) = 0
            END IF

            ! Deallocate arrays
            DEALLOCATE(temp_tn, temp_sn, temp_sshn)
         ELSE   
            ! Mark all as missing data
            tsn(:,:,:,:) = fbrmdi
            sshn(:,:) = fbrmdi 
         ENDIF
         ! Close netcdf file
         WRITE(numout,*) "Closing :", TRIM(filename)
         istat = nf90_close(ncid)
      END IF
   END SUBROUTINE ooo_read_file

   SUBROUTINE ooo_read_juld(filename, ifcst, julian)
      USE calendar
      !!--------------------------------------------------------------------
      !!                 *** ooo_read_juld ***
      !!
      !!   Purpose : To read model julian day information from file
      !!   Author  : A. Ryan Nov 2010
      !!--------------------------------------------------------------------

      !! Routine arguments
      CHARACTER(len=*), INTENT(IN)  :: filename
      INTEGER,          INTENT(IN)  :: ifcst
      REAL,             INTENT(OUT) :: julian    !: Julian day

      !! Local variables
      INTEGER :: year,  &   !: Date information
               & month, &
               & day,   &
               & hour,  &
               & minute,&
               & second
      INTEGER :: istat, &   !: Netcdf variables
               & ncid,  &
               & dimid, &
               & varid, &
               & ntime      
      REAL,DIMENSION(:),ALLOCATABLE :: r_sec     !: Remainder seconds
      CHARACTER(len=120) :: time_str  !: time string

      time_str=''
      !! Read in time_counter and remainder seconds
      istat = nf90_open(trim(filename),nf90_nowrite,ncid)
      istat = nf90_inq_dimid(ncid,'time_counter',dimid)
      IF (istat /= nf90_noerr) THEN
          istat = nf90_inq_dimid(ncid,'time',dimid)
      ENDIF
      istat = nf90_inquire_dimension(ncid,dimid,len=ntime)
      istat = nf90_inq_varid(ncid,'time_counter',varid)
      IF (istat /= nf90_noerr) THEN
          istat = nf90_inq_dimid(ncid,'time',dimid)
      ENDIF
      istat = nf90_get_att(ncid,varid,'units',time_str)
      ALLOCATE(r_sec(ntime))
      istat = nf90_get_var(ncid,varid, r_sec)
      istat = nf90_close(ncid)

      !! Fill yyyy-mm-dd hh:mm:ss
      !! format(('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2))
      100 format((14x, I4.4,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2))
      READ( time_str, 100 ) year, month, day, hour, minute, second

      CALL ymds2ju(year, month, day, r_sec(ifcst), julian)

      !! To take a comment from the ymds2ju subroutine

   !- In the case of the Gregorian calendar we have chosen to use
   !- the Lilian day numbers. This is the day counter which starts
   !- on the 15th October 1582.
   !- This is the day at which Pope Gregory XIII introduced the
   !- Gregorian calendar.
   !- Compared to the true Julian calendar, which starts some
   !- 7980 years ago, the Lilian days are smaler and are dealt with
   !- easily on 32 bit machines. With the true Julian days you can only
   !- the fraction of the day in the real part to a precision of
   !- a 1/4 of a day with 32 bits.
      
      !! The obs operator routines prefer to calculate Julian days since 
      !! 01/01/1950 00:00:00
      !! In order to convert to the 1950 version we must adjust by the number
      !! of days between 15th October 1582 and 1st Jan 1950

      julian = julian - 134123.
      
      DEALLOCATE(r_sec)
      
   END SUBROUTINE ooo_read_juld

END MODULE ooo_read