source: tags/nemo_dev_x3/NEMO/OPA_SRC/DTA/dtatem_fdir.h90 @ 105

   SUBROUTINE dtatem ( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dtatem  ***
      !!                     
      !! ** Purpose :   Reads temperature data (Levitus monthly data)
      !!
      !! ** Method  :   Read on unit numtdt the interpolated Levitus 
      !!      temperature onto the global grid.
      !!      Data begin at january.
      !!      The value is centered at the middle of month.
      !!      In the opa model, kt=1 agree with january 1.
      !!      At each time step, a linear interpolation is applied between
      !!      two monthly values.
      !!      Read on unit numtdt
      !!
      !! ** Action :
      !!      define t_dta array at time-step kt
      !!
      !! References :
      !!      Sydney Levitus, climatological atlas of the world ocean
      !!      NOAA professional paper 13, december 1982
      !!
      !! History :
      !!        !  91-03  ()  Original code
      !!        !  92-07  (M. Imbard)
      !!        !  99-10  (M.A. Foujols, M. Imbard)  NetCDF FORMAT
      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
      !!----------------------------------------------------------------------
      !! * Arguments
      INTEGER, INTENT( in ) ::   kt         ! ocean time-step

      !! * Local declarations
      INTEGER, PARAMETER ::   jpmois=12, jpf=1
      INTEGER ::   ji, jj, jl, ios, ik
      INTEGER ::   ibloc, ilindta
      INTEGER ::   imois, iman
      INTEGER ::   iimlu, ijmlu, ikmlu, ilmlu, immlu
      INTEGER ::   i15
      REAL(wp) ::   zxy, zl
      CHARACTER (len=20) ::   cltit
      !!----------------------------------------------------------------------


      ! 0. Initialization
      ! -----------------

      iman  = jpmois
      i15 = nday/16
      imois = nmonth + i15 - 1
      IF( imois == 0 ) imois = iman


      ! 1. First call kt=nit000
      ! -----------------------

      IF( kt == nit000 .AND. nlecte == 0 ) THEN
         ! open temp.dta file
         ibloc = 4096
         ilindta = ibloc*((jpidta*jpjdta*jpbytda-1 )/ibloc+1)
         CALL ctlopn(numtdt,'data_1m_potiential_temperature_nomask','OLD', 'UNFORMATTED', 'DIRECT',   &
                     ilindta,numout,lwp,1)
         ntem1 = 0
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'dta_tem : read monthly temperature in direct acces file'
         IF(lwp) WRITE(numout,*) '~~~~~~~'
         IF(lwp) WRITE(numout,*)

         ! 1.2 Read first records

         ! title, dimensions and tests

         READ( numtdt, REC=1, IOSTAT=ios ) cltit, iimlu, ijmlu, ikmlu,   &
                                           ilmlu, immlu
         IF( ios /= 0 ) THEN
            IF(lwp) WRITE(numout,*) 'e r r o r read numtdt ', ios
            STOP 'dtatem'
         ELSE
            IF ( iimlu /= jpidta ) STOP 4050
            IF ( ijmlu /= jpjdta ) STOP 4060
            IF ( ikmlu /= jpk    ) STOP 4070
            IF ( ilmlu /= jpmois ) STOP 4080
            IF ( immlu /= jpf    ) STOP 4090
            IF(lwp) WRITE(numout,*) 'nb of points in the 5 directions '
            IF(lwp) WRITE(numout,*) iimlu, ijmlu, ikmlu, ilmlu, immlu
         ENDIF

      ENDIF


      ! 2. Read monthly file
      ! -------------------

      IF( ( kt == nit000 .AND. nlecte == 0 ) .OR. imois /= ntem1 ) THEN
         nlecte = 1

         ! 2.1 Calendar computation
         
         ntem1 = imois        ! first file record used 
         ntem2 = ntem1 + 1    ! last  file record used
         ntem1 = MOD( ntem1, iman )
         IF( ntem1 == 0 )   ntem1 = iman
         ntem2 = MOD( ntem2, iman )
         IF( ntem2 == 0 )   ntem2 = iman
         IF(lwp) WRITE(numout,*) 'first record file used ntem1 ', ntem1
         IF(lwp) WRITE(numout,*) 'last  record file used ntem2 ', ntem2
         
         ! 2.3 Read monthly temperature data Levitus 

         CALL read3d( numtdt, temdta(1,1,1,1), ntem1+1 )
         CALL read3d( numtdt, temdta(1,1,1,2), ntem2+1 )
         
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) ' read Levitus temperature ok'
         IF(lwp) WRITE(numout,*)
         
         ! 2.4 Masks

         DO jl = 1, 2
            temdta(:,:,:,jl) = temdta(:,:,:,jl) * tmask(:,:,:)
            temdta(:,:,jpk,jl) = 0.
            IF( lk_zps ) THEN                ! z-coord. with partial steps
               DO jj = 1, jpj                     ! interpolation of temperature at the last level
                  DO ji = 1, jpi
                     ik = mbathy(ji,jj) - 1
                     IF ( ik > 2 ) THEN
                        zl = ( gdept(ik) - fsgdept(ji,jj,ik) ) / ( gdept(ik) - gdept(ik-1) )
                        temdta(ji,jj,ik,jl) = (1.-zl) * temdta(ji,jj,ik,jl) +zl * temdta(ji,jj,ik-1,jl
                     ENDIF
                  END DO
               END DO
            ENDIF
         END DO

         IF(lwp) THEN
            WRITE(numout,*) 'temperature Levitus month ', ntem1, ntem2
            WRITE(numout,*)
            WRITE(numout,*) ' Levitus mounth = ', ntem1, '  level = 1'
            CALL prihre( temdta(1,1,1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
            WRITE(numout,*) ' Levitus mounth = ', ntem1, '  level = ', jpk/2
            CALL prihre( temdta(1,1,jpk/2,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
            WRITE(numout,*) ' Levitus mounth = ',ntem1,'  level = ', jpkm1
            CALL prihre( temdta(1,1,jpkm1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
         ENDIF
      ENDIF

 
      ! 2. At every time step compute temperature data
      ! ----------------------------------------------

      zxy = FLOAT( nday + 15 - 30 * i15 ) / 30.
      t_dta(:,:,:) = (1.-zxy) * temdta(:,:,:,1) + zxy * temdta(:,:,:,2)

   END SUBROUTINE dta_tem