Index: /trunk/NEMO/OPA_SRC/SBC/flx_core.h90
===================================================================
--- /trunk/NEMO/OPA_SRC/SBC/flx_core.h90	(revision 472)
+++ /trunk/NEMO/OPA_SRC/SBC/flx_core.h90	(revision 472)
@@ -0,0 +1,956 @@
+   !!----------------------------------------------------------------------
+   !!                    ***  flx_core.h90  ***
+   !!----------------------------------------------------------------------
+   !!   flx     : update surface thermohaline fluxes from the NCAR data set
+   !!             read in a NetCDF file
+   !!----------------------------------------------------------------------
+   !! * Modules used     C A U T I O N  already defined in flxmod.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &
+      gsss,            &  !: SSS mean on nfbulk ocean time step
+      gu  , gv ,       &  !: (un,vn)(jk=1) mean over nfbulk ocean time step
+      !                   ! these variables are used for output in diawri
+      qlw_oce  ,       &  !: long wave flx over ocean
+      qla_oce  ,       &  !: latent heat flx over ocean
+      qsb_oce  ,       &  !: sensible heat flx over ocean
+      qlw_ice  ,       &  !: long wave flx over ice
+      qsb_ice             !: sensible heat flx over ice
+
+   !! * Module variables
+   INTEGER, PARAMETER ::   jpfile  = 8  ! maximum number of files to read 
+   CHARACTER (LEN=32), DIMENSION (jpfile) ::  clvarname  
+   CHARACTER (LEN=50), DIMENSION (jpfile) ::  clname 
+   INTEGER   ::      isnap
+   INTEGER, DIMENSION(jpfile) ::   &
+      numflxall,           &  ! logical units for surface fluxes data
+      nrecflx  , nrecflx2     ! first and second record to be read in flux file
+   REAL(wp), DIMENSION(jpfile) ::   &
+      freqh                 ! Frequency for each forcing file (hours)
+      !                     ! a negative value of -12 corresponds to monthly
+   REAL(wp), DIMENSION(jpi,jpj,jpfile,2) ::   &
+      flxdta                !:  set of NCAR 6hourly/daily/monthly fluxes
+   !!----------------------------------------------------------------------
+   !! History :
+   !!   9.0  !  04-08  (U. Schweckendiek)  Original code
+   !!        !  05-04  (L. Brodeau, A.M. Treguier)  severals modifications:
+   !!        !         - new bulk routine for efficiency
+   !!        !         - WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files
+   !!        !         - file names and file characteristics in namelist 
+   !!        !         - Implement reading of 6-hourly fields   
+   !!        !  06-02  (S. Masson, G. Madec)  IOM read + NEMO "style"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2006)
+   !! $Header$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE flx( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE flx  ***
+      !!                   
+      !! ** Purpose :   provide the thermohaline fluxes (heat and freshwater)
+      !!                and momentum fluxes (tau) 
+      !!      to the ocean at each time step.
+      !!
+      !! ** Method  :   Read NCAR data in a NetCDF file
+      !!               (file names and frequency of inputs specified in namelist)
+      !!          precipitation:  total (rain+snow)
+      !!          precipitation:  snow only 
+      !!          u10,v10 -> scalar wind at 10m in m/s -  ON 'T' GRID POINTS!!!
+      !!          solar radiation (short wave)  in W/m2
+      !!          thermal radiation (long wave) in W/m2
+      !!          specific humidity             in %
+      !!          temperature at 10m            in degrees K
+      !!
+      !!  ** Action  :   
+      !!          call flx_blk_albedo to compute ocean and ice albedo
+      !!          Calculates forcing fluxes to input into ice model 
+      !          or to be used directly for the ocean in ocesbc. 
+      !!
+      !!  ** Outputs 
+      !!          COMMENTS TO BE ADDED -units to be verified!
+      !!          taux    : zonal      wind stress on "u" points  (N/m2)
+      !!          tauy    : meridional wind stress on "v" points  (N/m2)
+      !!          qsr_oce : Solar flux over the ocean (W/m2)
+      !!          qnsr_oce: longwave flux over the ocean (W/m2)
+      !!          qsr_ice : solar flux over the ice (W/m2)
+      !!          qnsr_ice: longwave flux over the ice (W/m2)
+      !!          qla_ice : latent flux over sea-ice (W/m2)
+      !!          dqns_ice: total heat fluxes sensitivity over ice (dQ/dT)
+      !!          dqla_ice: latent heat flux sensitivity over ice (dQla/dT)
+      !!          fr1_i0
+      !!          fr2_i0
+      !!          tprecip
+      !!          sprecip
+      !!          evap
+      !!          
+      !! caution : now, in the opa global model, the net upward water flux is
+      !! -------   with mm/day unit,i.e. Kg/m2/s.
+      !!---------------------------------------------------------------------
+      !! * modules used
+      USE iom
+      USE par_oce
+      USE flx_oce
+      USE blk_oce         ! bulk variable
+      USE taumod
+      USE bulk
+      USE phycst
+      USE lbclnk
+      USE dtatem          ! FOR Ce = F(SST(levitus)):
+      
+      !! * arguments
+      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
+
+      !! * Local declarations
+      INTEGER , PARAMETER :: jpmaxtime = 366*4   ! maximum time steps in file will be for a 6 hourly field
+                                                 ! and a leap year (necessary variable for flinopen??)
+      INTEGER ::  irectot, irecflx
+      INTEGER ::  ihour                ! current hour of the day at which we read the forcings
+      INTEGER ::   &
+         imois, imois2, itime,      &  ! temporary integers
+         i15  , iman  ,             &  !    "          "
+         ifpr , jfpr  ,             &  ! frequency of index for print in prihre
+         jj   ,   ji                   !  Loop indices 
+      REAL(wp) ::   &
+         zadatrjb,                  &  ! date (noninteger day) at which we read the forcings
+         zsecond,                   &  ! temporary scalars
+         zxy    ,                   &  ! scalar for temporal interpolation
+         zcof                          ! scalar 
+      REAL(wp) ::   em , aa
+      REAL(wp) ::   zind1, zind2, zind3
+      REAL(wp), DIMENSION(jpi,jpj, jpfile) ::   &
+         flxnow             ! input flux values at current time step 
+      REAL(wp), DIMENSION(jpi,jpj) ::   &
+         sstk  ,   t04  
+      REAL(wp), DIMENSION(jpi,jpj) ::   &
+         qtune    ,      &  ! artifical field for tuning q
+         dqlw_ice ,      &  ! long wave heat flx sensitivity over ice
+         dqsb_ice ,      &  ! sensible heat flx sensitivity over ice
+         alb_ice  ,      &  ! albedo of ice
+         alb_oce_os,     &  ! albedo of the ocean under overcast sky
+         alb_ice_os,     &  ! albedo of the ice under overcast sky
+         alb_ice_cs,     &  ! albedo of ice under clear sky
+         alb_oce_cs,     &  ! albedo of the ocean under clear sky
+         zsstlev,        &  ! SST from Levitus in K
+         zsst,           &  ! nfbulk : mean SST
+         zsss,           &  ! nfbulk : mean tn_ice(:,:,1)
+         zut,            &  ! nfbulk : mean U at T-point
+         zvt,            &  ! nfbulk : mean V at T-point
+         dUnormt,        &  ! scalar wind (norm) on T points
+         tauxt, tauyt,   &  ! wind stress computed at T-point
+         qsatw,          &  ! specific humidity at zSST
+         qsat,           &  ! specific humidity at zsss 
+         Ch,             &  ! coefficient for sensible heat flux
+         Ce,             &  ! coefficient for latent heat flux
+         Cd                 ! coefficient for wind stress
+      REAL(wp), PARAMETER ::   &
+         rhoa  = 1.22,    &  ! air density
+         cp    = 1000.5,  &  ! specific heat of air
+         Lv    = 2.5e6,   &  ! latent heat of vaporization
+         Ls    = 2.839e6, &  ! latent heat of sublimation
+         Stef  = 5.67e-8, &  ! Stefan Boltzmann constant
+         Cice = 1.63e-3      ! transfert coefficient over ice
+      REAL(wp), DIMENSION(jpi,jpj) ::   &
+         catm1               ! 
+
+      NAMELIST /namcore/  clname,  clvarname, freqh
+      !!---------------------------------------------------------------------
+
+      !! =============================================== !!
+      !!   PART A:   READING FLUX FILES WHEN NECESSARY   !!
+      !! =============================================== !!
+
+      !--- calculation for monthly data 
+      i15 = INT( 2 * FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) )
+      iman  = 12
+      imois = nmonth + i15 - 1
+      IF( imois == 0 ) imois = iman
+      imois2 = nmonth
+
+      !--- calculation for 6-hourly data 
+      !    we need the fraction of day, measured in hours at the date of forcings. 
+      !    This is the time step before the date we are calculating, zadatrj 
+      !    "adatrj" (real time in days (noninteger))
+      zadatrjb = adatrj - rdt / rday
+      ihour = INT( ( zadatrjb - INT(zadatrjb) ) * 24 )
+      
+      !                                         ! ------------------------ !
+      IF( kt == nit000 ) THEN                   !   first call kt=nit000   !
+         !                                      ! ------------------------ !
+
+         !--- Initializes default values of file names, frequency of forcings 
+         !    and variable to be read in the files, before reading namelist 
+         clname(1) = 'precip_core.nc'   ; freqh(1) = -12 ; clvarname(1) = 'precip'     ! monthly
+         clname(2) = 'u10_core.nc'      ; freqh(2) =  24 ; clvarname(2) = 'u_10'       ! 6 hourly          
+         clname(3) = 'v10_core.nc'      ; freqh(3) =  24 ; clvarname(3) = 'v_10'       ! 6 hourly    
+         clname(4) = 'q10_core.nc'      ; freqh(4) =  24 ; clvarname(4) = 'q_10'       ! 6 hourly   
+         clname(5) = 'tot_solar_core.nc'; freqh(5) =  24 ; clvarname(5) = 'tot_solar'  ! daily 
+         clname(6) = 'therm_rad_core.nc'; freqh(6) =  24 ; clvarname(6) = 'therm_rad'  ! daily 
+         clname(7) = 'temp_10m_core.nc' ; freqh(7) =  24 ; clvarname(7) = 't_10'       ! 6 hourly 
+         clname(8) = 'snow_core.nc'     ; freqh(8) = -12 ; clvarname(8) = 'snow'       ! monthly
+
+         !--- Read Namelist namcore : OMIP/CORE 
+         REWIND ( numnam )
+         READ   ( numnam, namcore )
+         IF(lwp) THEN
+            WRITE(numout,*)' '
+            WRITE(numout,*)' flxmod/flx_core : global CORE fields in NetCDF format'
+            WRITE(numout,*)' ~~~~~~~~~~~~~~~ '
+            WRITE(numout,*) '      list of files and frequency (hour), or monthly (-12) '
+            DO ji = 1, jpfile
+                WRITE(numout,*) trim(clname(ji)), ' frequency:', freqh(ji)
+            END DO
+         ENDIF
+
+         !--- Initialization to zero
+         qsr_oce (:,:) = 0.e0
+         qsb_oce (:,:) = 0.e0
+         qla_oce (:,:) = 0.e0
+         qlw_oce (:,:) = 0.e0
+         qnsr_oce(:,:) = 0.e0
+         qsr_ice (:,:) = 0.e0
+         qnsr_ice(:,:) = 0.e0
+         qla_ice (:,:) = 0.e0
+         qlw_ice (:,:) = 0.e0
+         qsb_ice (:,:) = 0.e0
+
+         dqns_ice(:,:) = 0.e0
+         dqla_ice(:,:) = 0.e0
+         tprecip (:,:) = 0.e0
+         sprecip (:,:) = 0.e0
+         evap    (:,:) = 0.e0
+
+         flxnow  (:,:,:) = 1.e-15         !RB
+         flxdta  (:,:,:,:) = 1.e-15       !RB
+
+         nrecflx (:)    = 0      ! switch for reading flux data for each file
+         nrecflx2(:)    = 0      ! switch for reading flux data for each file
+         
+         !--- Open the files of the list
+         DO ji = 1, jpfile
+            CALL iom_open( clname(ji), numflxall(ji) )
+         END DO
+
+      ENDIF
+
+      !                                         ! ------------------------ !
+      !                                         !  Read files if required  !
+      !                                         ! ------------------------ !
+
+      !--- read data if necessary checks each file in turn 
+
+      DO ji = 1, jpfile
+
+         IF (  freqh(ji) .GT.0 .AND.  freqh(ji) .LT. 24 ) THEN          !--- Case of 6-hourly flux data  
+            !--- calculate current snapshot from hour of day 
+            isnap = ihour / INT( freqh(ji) ) + 1
+            !--- reading is necessary when nrecflx(ji) differs from isnap 
+            IF( nrecflx(ji) /= isnap ) THEN   
+               nrecflx(ji) = isnap
+               irecflx = (nday_year-1)*24 / freqh(ji) + isnap
+               irectot = 365*24 / freqh(ji)        !  this variable is not used by iom_get
+               IF(lwp) WRITE (numout,*)' CORE flux 6-h record :',irecflx, ' var:',trim(clvarname(ji))
+               CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), irecflx )
+            ENDIF
+
+         ELSE IF ( freqh(ji) .EQ. 24 ) THEN                             !--- Case of daily flux data  
+
+            !--- reading is necessary when nrecflx(ji) differs from nday
+            IF( nrecflx(ji) /= nday ) THEN   
+               nrecflx(ji) = nday          !! remember present read day
+               irecflx     = nday_year     !! 
+               irectot = 365        !  this variable is not used by iom_get
+               IF(lwp) WRITE (numout,*)'DAILY CORE flux record :',irecflx, ' var:',trim(clvarname(ji))
+               CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), irecflx )
+            ENDIF
+
+         ELSE IF ( freqh(ji) .EQ. -12 ) THEN                            !--- Case monthly data from CORE
+            ! we read two months all the time although we 
+            ! could only read one and swap the arrays 
+            IF( kt == nit000 .OR. imois /= nrecflx(ji) ) THEN
+               ! calendar computation
+               ! nrecflx number of the first file record used in the simulation
+               ! nrecflx2 number of the last  file record
+               nrecflx(ji) = imois
+               nrecflx2(ji) = nrecflx(ji)+1
+               nrecflx(ji) = MOD( nrecflx(ji), iman )
+               IF( nrecflx(ji) == 0 )   nrecflx(ji) = iman
+               nrecflx2(ji) = MOD( nrecflx2(ji), iman )
+               IF ( nrecflx2(ji) == 0 ) nrecflx2(ji) = iman
+               irectot = 12        !  this variable is not used by iom_get
+               IF(lwp) WRITE(numout,*) 'MONTHLY CORE flux record 1:', nrecflx(ji),   &
+                  &                    ' rec 2:', nrecflx2(ji), ' var:', trim(clvarname(ji))
+               CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), nrecflx (ji) )
+               CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,2), nrecflx2(ji) )
+            ENDIF
+         ENDIF
+      END DO    !   end of loop over forcing files to verify if reading is necessary 
+ 
+      !--- Printout if required      
+      !  IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN
+      !            printout at the initial time step only (unless you want to debug!!)
+      IF( kt == nit000 ) THEN
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) ' read daily and monthly CORE fluxes: ok'
+            WRITE(numout,*)
+            ifpr = INT(jpi/8)      ;      jfpr = INT(jpj/10)
+            DO ji = 1, jpfile 
+               WRITE(numout,*) trim(clvarname(ji)),' day: ',ndastp
+               CALL prihre(flxdta(1,1,ji,1),jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout)
+               WRITE(numout,*)
+            END DO
+         ENDIF
+      ENDIF
+
+      !                                         ! ------------------------ !
+      !                                         !   Interpolates in time   !
+      !                                         ! ------------------------ !
+      ! N.B. For now, only monthly values are interpolated, and they 
+      !      are interpolated to the current day, not to the time step.
+      !
+      DO ji = 1, jpfile       
+         IF( freqh(ji) .EQ. -12 ) THEN
+            zxy = REAL(nday) / REAL( nobis(imois) ) + 0.5 - i15         !!! <== Caution nobis hard coded !!!
+            flxnow(:,:,ji) = ( ( 1. - zxy) * flxdta(:,:,ji,1) + zxy * flxdta(:,:,ji,2) )
+         ELSE
+            flxnow(:,:,ji) = flxdta(:,:,ji,1)
+         ENDIF
+      END DO
+      ! JMM : add vatm needed in tracer routines  ==> wind module ???
+      vatm(:,:) = SQRT( flxnow(:,:,2) * flxnow(:,:,2) + flxnow(:,:,3) * flxnow(:,:,3) )
+
+
+      !! =============================================== !!
+      !!   PART B:   CORE BULK CALCULATION               !!
+      !! =============================================== !!
+
+      !  for other forcing cases this is done in modules bulk.F90 and flxblk
+
+      !                                         ! ------------------------ !
+      !                                         !   Bulk initialisation    !
+      !                                         ! ------------------------ !
+
+      ! code part from bulk.F90 :
+
+      IF( kt == nit000) THEN
+         !--- computation of rdtbs2     ===>  !gm  is it really usefull ????
+         IF( nacc == 1 ) THEN
+            rdtbs2 = nfbulk * rdtmin * 0.5
+         ELSE
+           rdtbs2 = nfbulk * rdt * 0.5
+         ENDIF
+         IF( .NOT.ln_rstart ) THEN
+            zcof =  REAL( nfbulk - 1, wp )
+            gsst(:,:) =  zcof * ( tn(:,:,1)  + rt0 )
+            gsss(:,:) =  zcof * tn_ice(:,:)
+            gu  (:,:) =  zcof * un(:,:,1)
+            gv  (:,:) =  zcof * vn(:,:,1)
+         ENDIF
+      ENDIF
+
+      ! ----------------------------------------------------------------------------- !
+      !      0   Mean SST, ice temperature and ocean currents over nfbulk pdt         !
+      ! ----------------------------------------------------------------------------- !
+      !
+      gsst(:,:) = gsst(:,:) + (tn(:,:,1) + rt0 )  
+      gsss(:,:) = gsss(:,:) + tn_ice(:,:)
+      gu  (:,:) = gu  (:,:) + un    (:,:,1)
+      gv  (:,:) = gv  (:,:) + vn    (:,:,1)
+
+      IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN
+         !
+         zsst(:,:) = gsst(:,:) / REAL( nfbulk, wp ) * tmask(:,:,1) ! mean sst in K
+         zsss(:,:) = gsss(:,:) / REAL( nfbulk ) * tmask(:,:,1) ! mean tn_ice in K
+
+         where ( zsst(:,:) .eq. 0.e0 ) zsst(:,:) = rt0     !lb vilain !!???
+         where ( zsss(:,:) .eq. 0.e0 ) zsss(:,:) = rt0     !lb    //
+
+!!gm  better coded:
+         ! Caution set to rt0 over land, not 0 Kelvin (otherwise floating point exception in bulk computation) 
+!        zcof =  1. / REAL( nfbulk, wp )
+!        zsst(:,:) = gsst(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1)    ! mean sst in K
+!        zsss(:,:) = gsss(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1)   ! mean tn_ice in K
+!!gm
+       
+         zut(:,:) = 0.e0      !!gm not necessary but at least first and last column
+         zvt(:,:) = 0.e0
+         !            lb
+         ! Interpolation of surface current at T-point, zut and zvt :
+         DO ji=2, jpi
+            zut(ji,:) = 0.5*(gu(ji-1,:) + gu(ji,:)) / REAL( nfbulk )
+         END DO
+         !
+         DO jj=2, jpj
+            zvt(:,jj) = 0.5*(gv(:,jj-1) + gv(:,jj)) / REAL( nfbulk )
+         END DO
+!!gm better coded
+!        zcof =  0.5 / REAL( nfbulk, wp )
+!        DO jj = 2, jpjm1
+!           DO ji = fs_2, fs_jpim1   ! vector opt.
+!           zut(ji,jj) = ( gu(ji-1,jj  ) + gu(ji,jj) ) * zcof
+!           zvt(ji,jj) = ( gv(ji  ,jj-1) + gv(ji,jj) ) * zcof
+!           END DO
+!        END DO
+!!gm
+
+         CALL lbc_lnk( zut(:,:), 'T', -1. )         
+         CALL lbc_lnk( zvt(:,:), 'T', -1. )         
+
+         ! ----------------------------------------------------------------------------- !
+         !      I   Radiative FLUXES                                                     !
+         ! ----------------------------------------------------------------------------- !
+
+         !--- Ocean & Ice Albedo
+         alb_oce_os(:,:) = 0.   ;   alb_oce_cs(:,:) = 0.   !gm is it necessary ???
+         alb_ice_os(:,:) = 0.   ;   alb_ice_cs(:,:) = 0.
+         CALL flx_blk_albedo( alb_ice_os, alb_oce_os, alb_ice_cs, alb_oce_cs )
+
+         !--- Radiative fluxes over ocean
+         qsr_oce(:,:) = ( 1. - 0.066 ) * flxnow(:,:,5)                                 ! Solar Radiation
+         qlw_oce(:,:) = flxnow(:,:,6) - Stef*zsst(:,:)*zsst(:,:)*zsst(:,:)*zsst(:,:)   ! Long Waves (Infra-red)
+
+         !--- Radiative fluxes over ice
+         qsr_ice(:,:) = ( 1. - alb_ice_cs(:,:) ) * flxnow(:,:,5)                                  ! Solar Radiation
+         qlw_ice(:,:) = 0.95 * ( flxnow(:,:,6) - Stef*zsss(:,:)*zsss(:,:)*zsss(:,:)*zsss(:,:) )   ! Long Waves (Infra-red)
+         !-------------------------------------------------------------------------------!
+       
+       
+         ! ----------------------------------------------------------------------------- !
+         !     II    Turbulent FLUXES                                                    !
+         ! ----------------------------------------------------------------------------- !
+         !
+         ! scalar wind ( = | U10m - SSU | )
+         ! It is important to take into account the sea surface courant
+         !                lb
+         ! Now, wind components are provided on T-points within the netcdf input file.
+         ! Thus, the wind module is computded at T-points taking into account the sea
+
+         !--- Module of relative wind
+         dUnormt(:,:) = sqrt( (flxnow(:,:,2) - zut(:,:))*(flxnow(:,:,2) - zut(:,:))  &
+              &             + (flxnow(:,:,3) - zvt(:,:))*(flxnow(:,:,3) - zvt(:,:)) )
+!!gm  more efficient coding:
+!        DO jj = 1, jpi
+!           DO ji = 1, jpi
+!              zzu = flxnow(ji,jj,2) - zut(ji,jj)
+!              zzv = flxnow(ji,jj,3) - zvt(ji,jj)
+!              dUnormt(ji,jj) = SQRT( zzu*zzu + zzv*zzv )
+!           END DO 
+!        END DO 
+!!gm
+         !                lb
+         !
+         !--- specific humidity at temp SST
+         qsatw(:,:) = 0.98*640380*exp(-5107.4/zsst(:,:))/rhoa
+         !
+         !--- specific humidity at temp tn_ice
+         qsat(:,:)  = 11637800*exp(-5897.8/zsss(:,:))/rhoa
+         !
+         !--- NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point (NCAR Bulk formulae)
+         CALL TURB_CORE( 10., zsst(:,:), flxnow(:,:,7), qsatw(:,:), flxnow(:,:,4),   &
+            &                              dUnormt(:,:), Cd(:,:), Ch(:,:), Ce(:,:) )
+  
+         !   II.1) Momentum over ocean and ice
+         !   ---------------------------------
+         ! Tau_x at T-point
+         tauxt(:,:) = rhoa*dUnormt(:,:)*( (1. - freeze(:,:))*Cd(:,:)*(flxnow(:,:,2) - zut(:,:))    &
+              &                          + freeze(:,:)*Cice*flxnow(:,:,2) ) !lb correct pour glace
+         ! Tau_y at T-point
+         tauyt(:,:) = rhoa*dUnormt(:,:)*( (1. - freeze(:,:))*Cd(:,:)*(flxnow(:,:,3) - zvt(:,:))    &
+              &                          + freeze(:,:)*Cice*flxnow(:,:,3) )
+         !
+         CALL lbc_lnk( tauxt(:,:), 'T', -1. )       !!gm seems not decessary at this point....
+         CALL lbc_lnk( tauyt(:,:), 'T', -1. )
+         !
+         !Tau_x at U-point
+         DO ji = 1, jpim1
+            taux(ji,:) = 0.5*(tauxt(ji,:) + tauxt(ji+1,:))
+         END DO
+         !
+         ! Tau_y at V-point
+         DO jj = 1, jpjm1
+            tauy(:,jj) = 0.5*(tauyt(:,jj) + tauyt(:,jj+1))
+         END DO
+
+!!gm ==>  at this point, a lbc is required, no?
+         !
+         ! lb : should we do this here?    !!gm yes should we remove that???
+         tauxg(:,:) = taux(:,:)  ! Save components in 
+         tauyg(:,:) = tauy(:,:)  ! geographical ref on U grid
+         !
+         !
+         !   II.2) Turbulent fluxes over ocean
+         !   ---------------------------------
+         !
+         ! Sensible Heat :
+         qsb_oce(:,:) = rhoa * cp * Ch(:,:) * ( zsst(:,:) - flxnow(:,:,7) ) * dUnormt(:,:)
+         !
+         ! Latent Heat :
+         evap(:,:)    = rhoa * Ce(:,:) * ( qsatw(:,:) - flxnow(:,:,4) ) * dUnormt(:,:)
+         qla_oce(:,:) = Lv * evap(:,:)
+         !
+         ! 
+         !   II.3) Turbulent fluxes over ice
+         !   -------------------------------
+         ! 
+         ! Sensible Heat :
+         qsb_ice(:,:) = rhoa*cp*Cice*(zsss(:,:) - flxnow(:,:,7))*dUnormt(:,:) !lb     use
+         !                                                                    !lb   dUnormt???
+         ! Latent Heat :                                                      !lb or rather Unormt?
+         qla_ice(:,:) = Ls*rhoa*Cice*(qsat(:,:) - flxnow(:,:,4))*dUnormt(:,:)
+         !
+         !-------------------------------------------------------------------------------------
+         
+         
+         
+         ! ----------------------------------------------------------------------------- !
+         !     III    Total FLUXES                                                       !
+         ! ----------------------------------------------------------------------------- !
+         !
+         !   III.1) Downward Non Solar flux over ocean
+         !   -----------------------------------------
+         qnsr_oce(:,:) = qlw_oce(:,:) - qsb_oce(:,:) - qla_oce(:,:)
+         !
+         !   III.1) Downward Non Solar flux over ice
+         !   ---------------------------------------
+         qnsr_ice(:,:) = qlw_ice(:,:) - qsb_ice(:,:) - qla_ice(:,:) 
+         !
+         !-------------------------------------------------------------------------------!
+         
+         
+         
+         ! 6. TOTAL NON SOLAR SENSITIVITY 
+       
+         dqlw_ice(:,:)= 4.0*0.95*Stef*zsss(:,:)*zsss(:,:)*zsss(:,:)
+         
+         ! d qla_ice/ d zsss
+         dqla_ice(:,:) = -Ls*Cice*0.98*11637800/(rhoa*rhoa) &
+              * (-5897.8)/(zsss(:,:)*zsss(:,:)) &
+              * exp(-5897.8/zsss(:,:)) * dUnormt(:,:)
+         
+         ! d qsb_ice/ d zsss
+         dqsb_ice(:,:) = rhoa * cp * Cice * dUnormt(:,:)
+         
+         dqns_ice(:,:) = - ( dqlw_ice(:,:) + dqsb_ice(:,:) + dqla_ice(:,:) )
+         
+         
+         !--------------------------------------------------------------------
+         ! FRACTION of net shortwave radiation which is not absorbed in the
+         ! thin surface layer and penetrates inside the ice cover
+         ! ( Maykut and Untersteiner, 1971 ; Elbert and Curry, 1993 )
+         
+         catm1(:,:) = 1.0  - 0.3  ! flxnow(:,:,8)
+         fr1_i0(:,:) = & 
+              (0.18 * catm1(:,:) + 0.35 * flxnow(:,:,8))
+         fr2_i0(:,:) = &
+              (0.82 * catm1(:,:) + 0.65 * flxnow(:,:,8))
+         
+         
+         !  Total  PRECIPITATION (kg/m2/s)
+         tprecip(:,:) = flxnow(:,:,1) 
+         ! rename precipitation for freshwater budget calculations:
+         watm(:,:) =  flxnow(:,:,1)*1000
+         !
+         
+         !  SNOW PRECIPITATION  (kg/m2/s)
+         sprecip(:,:) =  flxnow(:,:,8) 
+       
+         !---------------------------------------------------------------------
+         
+         
+         
+         CALL lbc_lnk( taux (:,:)    , 'U', -1. )
+         CALL lbc_lnk( tauy (:,:)    , 'V', -1. )
+         CALL lbc_lnk( qsr_oce (:,:) , 'T', 1. )
+         CALL lbc_lnk( qnsr_oce(:,:) , 'T', 1. )
+         CALL lbc_lnk( qsr_ice (:,:) , 'T', 1. )
+         CALL lbc_lnk( qnsr_ice(:,:) , 'T', 1. )
+         CALL lbc_lnk( qla_ice (:,:) , 'T', 1. )
+         CALL lbc_lnk( dqns_ice(:,:) , 'T', 1. )
+         CALL lbc_lnk( dqla_ice(:,:) , 'T', 1. )
+         CALL lbc_lnk( fr1_i0  (:,:) , 'T', 1. )
+         CALL lbc_lnk( fr2_i0  (:,:) , 'T', 1. )
+         CALL lbc_lnk( tprecip (:,:) , 'T', 1. )
+         CALL lbc_lnk( sprecip (:,:) , 'T', 1. )
+         CALL lbc_lnk( evap    (:,:) , 'T', 1. )
+         !!
+         !!
+         !!   NEVER mask the windstress!!
+         qsr_oce (:,:) = qsr_oce (:,:)*tmask(:,:,1)
+         qnsr_oce(:,:) = qnsr_oce(:,:)*tmask(:,:,1)
+         qsr_ice (:,:) = qsr_ice (:,:)*tmask(:,:,1)
+         qnsr_ice(:,:) = qnsr_ice(:,:)*tmask(:,:,1)
+         qla_ice (:,:) = qla_ice (:,:)*tmask(:,:,1)
+         dqns_ice(:,:) = dqns_ice(:,:)*tmask(:,:,1)
+         dqla_ice(:,:) = dqla_ice(:,:)*tmask(:,:,1)
+         fr1_i0  (:,:) = fr1_i0  (:,:)*tmask(:,:,1)
+         fr2_i0  (:,:) = fr2_i0  (:,:)*tmask(:,:,1)
+         tprecip (:,:) = tprecip (:,:)*tmask(:,:,1)
+         sprecip (:,:) = sprecip (:,:)*tmask(:,:,1)
+         evap    (:,:) = evap    (:,:)*tmask(:,:,1)
+         gsst(:,:) = 0.e0
+         gsss(:,:) = 0.e0
+         gu  (:,:) = 0.e0
+         gv  (:,:) = 0.e0
+  
+         !   Degugging print (A.M. Treguier)       
+         !      write (55) taux, tauy, qnsr_oce, qsb_ice, qnsr_ice, qla_ice, dqns_ice &
+         !  &     , dqla_ice, fr1_i0, fr2_i0, qlw_oce, qla_oce, qsb_oce, evap
+         !      write(numout,*) 'written 14 arrays on unit fort.55'
+  
+         
+      ENDIF
+    
+      ! ------------------- !
+      ! Last call kt=nitend !
+      ! ------------------- !
+    
+      ! Closing of the numflx file (required in mpp)
+    
+      IF( kt == nitend ) THEN
+         DO ji=1, jpfile 
+            CALL iom_close(numflxall(ji))
+         END DO
+      ENDIF
+    
+   END SUBROUTINE flx
+  
+  
+   SUBROUTINE TURB_CORE( zzu, T_0, T_a, q_sat, q_a,   &
+      &                       dU , C_d, C_h  , C_e  )
+      !!----------------------------------------------------------------------
+      !!                      ***  ROUTINE  turb_core  ***
+      !!
+      !! ** Purpose :   Computes turbulent transfert coefficients of surface 
+      !!      fluxes according to Large & Yeager (2004)
+      !!
+      !! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D
+      !!      Momentum, Latent and sensible heat exchange coefficients
+      !!      Caution: this procedure should only be used in cases when air
+      !!      temperature (T_air), air specific humidity (q_air) and wind (dU)
+      !!      are provided at the same height 'zzu'!
+      !!
+      !! References :
+      !!      Large & Yeager, 2004 : ???
+      !! History :
+      !!        !  XX-XX  (???     )  Original code
+      !!   9.0  !  05-08  (L. Brodeau)  Optimisation
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      REAL(wp), INTENT( in  ) ::   &
+         zzu              ! height for all given atmospheric variables  [m]
+      REAL(wp), INTENT( in  ), DIMENSION(jpi,jpj) ::  &
+         T_0,       &     ! sea surface temperature                     [Kelvin]
+         T_a,       &     ! potential air temperature at zzu            [Kelvin]
+         q_sat,     &     ! sea surface specific humidity at zzu        [kg/kg]
+         q_a,       &     ! specific air humidity at zzu                [kg/kg]
+         dU               ! relative wind module |U(zzu)-U(0)| at zzu   [m/s]
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)  :: &
+         C_d,    &        ! transfer coefficient for momentum       (tau)
+         C_h,    &        ! transfer coefficient for sensible heat  (Q_sens)
+         C_e              ! tansfert coefficient for evaporation    (Q_lat)
+
+      !! * Local declarations
+      REAL(wp), DIMENSION(jpi,jpj)  ::   &
+         dU10,        &       ! dU                                   [m/s]
+         dT,          &       ! air/sea temperature differeence      [K]
+         dq,          &       ! air/sea humidity difference          [K]
+         Cd_n10,      &       ! 10m neutral drag coefficient
+         Ce_n10,      &       ! 10m neutral latent coefficient
+         Ch_n10,      &       ! 10m neutral sensible coefficient
+         Cd,          &       ! drag coefficient
+         Ce,          &       ! latent coefficient
+         Ch,          &       ! sensible coefficient
+         sqrt_Cd_n10, &       ! root square of Cd_n10
+         sqrt_Cd,     &       ! root square of Cd
+         T_vpot,      &       ! virtual potential temperature        [K]
+         T_star,      &       ! turbulent scale of tem. fluct.
+         q_star,      &       ! turbulent humidity of temp. fluct.
+         U_star,      &       ! turb. scale of velocity fluct.
+         L,           &       ! Monin-Obukov length                  [m]
+         zeta,        &       ! stability parameter at height zzu
+         X2, X,       &   
+         psi_m,       &
+         psi_h,       &
+         U_n10,       &       ! neutral wind velocity at 10m        [m]   
+         xlogt
+ 
+      INTEGER                    :: jkk,ji,jj,ii,ij,ilocu(2),  jk          ! doomy loop for iterations
+      INTEGER, PARAMETER         :: nit = 3       ! number of iterations
+      REAL(wp), DIMENSION(jpi,jpj) :: dbg1,dbg2,dbg3,dbg4
+      REAL :: zXmax
+
+      INTEGER, DIMENSION(jpi,jpj)  ::   &
+         &                  stab,        &       ! stability test integer
+         &                  stabit               ! stability within iterative loop
+
+      REAL(wp), PARAMETER ::                        &
+         &                     pi     = 3.14159,      &  ! Pi
+         &                     grav   = 9.8,          &  ! gravity
+         &                     kappa  = 0.4              ! von Karman's constant
+      !!----------------------------------------------------------------------
+      !!
+      !!                  -------------
+      !!                    S T A R T
+      !!                  -------------
+      !! 
+      !! I.1 ) Preliminary stuffs
+      !! ------------------------
+      !!
+      !! Air/sea differences
+      !! -------------------
+      dU10 = max(0.5, dU)     ! we don't want to fall under 0.5 m/s
+      dT   = T_a - T_0        ! assuming that T_a is allready the potential temp. at zzu
+      dq   = q_a - q_sat
+      !!    
+      !! Virtual potential temperature
+      !! -----------------------------
+      T_vpot = T_a*(1. + 0.608*q_a)
+      !!
+      !!
+      !! I.2 ) Computing Neutral Drag Coefficient
+      !! ----------------------------------------
+      Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !  \\ L & Y eq. (6a)
+      sqrt_Cd_n10 = sqrt(Cd_n10)
+      !!
+    Ce_n10  = 1E-3 * ( 34.6 * sqrt_Cd_n10 )               !  \\ L & Y eq. (6b)
+    !!
+    !! First guess of stabilitty :
+    stab    = 0.5 + sign(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0 
+    Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) !  \\ L & Y eq. (6c), (6d)
+    !!
+    !! Initializing transfert coefficients with their first guess neutral equivalents :
+    Cd = Cd_n10  ;  Ce = Ce_n10  ;  Ch = Ch_n10
+    !!
+    !!
+    !!
+    !! II ) Now starting iteration loop (IDM)
+    !! ----------------------------------------
+    !!
+    DO jk=1, nit
+       !!
+       sqrt_Cd = sqrt(Cd)
+       !!
+       !! Turbulent scales :
+       !! ------------------
+       U_star  = sqrt_Cd*dU10                !  \\ L & Y eq. (7a)
+       T_star  = Ch/sqrt_Cd*dT               !  \\ L & Y eq. (7b)
+       q_star  = Ce/sqrt_Cd*dq               !  \\ L & Y eq. (7c)
+       !!
+       !! Estimate the Monin-Obukov length :
+       !! ----------------------------------
+       !dbg1 = T_star/T_vpot 
+       !dbg2 = q_star/(q_a + 1./0.608) 
+       L  = (U_star**2)/( kappa*grav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) )
+       !!
+       !! Stability parameters :
+       !! ----------------------
+       zeta = zzu/L
+       zeta   = sign( min(abs(zeta),10.0), zeta )
+       !!
+       !! Psis, L & Y eq. (8c), (8d), (8e) :
+       !! ----------------------------------
+       X2 = sqrt(abs(1. - 16.*zeta))  ;  X2 = max(X2 , 1.0) ;  X  = sqrt(X2)
+       !!
+       stabit    = 0.5 + sign(0.5,zeta)
+       !!
+!      dbg1 =  -5*zeta*stabit
+!      dbg2 =  2*log((1. + X)/2)
+!      dbg3 = log((1. + X2)/2)
+!      dbg4 =  2*atan(X)
+
+       psi_m = -5*zeta*stabit  &                                                       ! Stable
+            & + (1 - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable
+       !!
+       psi_h = -5*zeta*stabit  &                                                       ! Stable
+            & + (1 - stabit)*(2*log( (1. + X2)/2 ))                                    ! Unstable
+       !!
+       !!
+       !! Shifting the wind speed to 10m and neutral stability :
+       !! ------------------------------------------------------
+       U_n10 = dU10*1./(1. + sqrt(Cd_n10)/kappa*(log(zzu/10.) - psi_m))    ! \\ L & Y eq. (9a)
+       !!
+       !!
+       !! Updating the neutral 10m transfer coefficients :
+       !! ------------------------------------------------
+       Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)        ! \\ L & Y eq. (6a)
+       sqrt_Cd_n10 = sqrt(Cd_n10)
+       !!
+       Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                     ! \\ L & Y eq. (6b)
+       !!
+       stab    = 0.5 + sign(0.5,zeta)
+       Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab))       ! \\ L & Y eq. (6c), (6d)
+       !!
+       !!
+       !! Shifting the neutral  10m transfer coefficients to ( zzu , zeta ) :
+       !! --------------------------------------------------------------------
+       !! Problem here, formulation used within L & Y differs from the one provided 
+       !! in their fortran code (only for Ce and Ch)
+       !!
+       Cd = Cd_n10/(1. + sqrt_Cd_n10/kappa*(log(zzu/10) - psi_m))**2     ! \\ L & Y eq. (10a)
+       !!
+       xlogt = log(zzu/10) - psi_h
+       !?      Ch = Ch_n10*sqrt(Cd/Cd_n10)/(1. + Ch_n10/(kappa*sqrt_Cd_n10)*xlogt)
+       Ch = Ch_n10/( 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 )**2             ! \\ L & Y eq. (10b)
+       !!
+       !?      Ce = Ce_n10*sqrt(Cd/Cd_n10)/(1. + Ce_n10/(kappa*sqrt_Cd_n10)*xlogt)
+       Ce = Ce_n10/( 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10 )**2             ! \\ L & Y eq. (10c)
+       !!
+       !!
+    END DO
+    !!
+    C_d = Cd
+    C_h = Ch
+    C_e = Ce
+    !!
+  END SUBROUTINE TURB_CORE
+  
+    
+  SUBROUTINE flx_blk_albedo( palb , palcn , palbp , palcnp )
+    !!----------------------------------------------------------------------
+    !!               ***  ROUTINE flx_blk_albedo  ***
+    !!
+    !! ** Purpose :   Computation of the albedo of the snow/ice system
+    !!      as well as the ocean one
+    !!
+    !! ** Method  : - Computation of the albedo of snow or ice (choose the
+    !!      right one by a large number of tests
+    !!              - Computation of the albedo of the ocean
+    !!
+    !! References :
+    !!      Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250.
+    !!
+    !! History :
+    !!  8.0   !  01-04  (LIM 1.0)
+    !!  8.5   !  03-07  (C. Ethe, G. Madec)  Optimization (old name:shine)
+    !!----------------------------------------------------------------------
+    !! * Modules used
+    USE ice                   ! ???
+    
+    !! * Arguments
+    REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::  &
+         palb         , & !  albedo of ice under overcast sky
+         palcn        , & !  albedo of ocean under overcast sky
+         palbp        , & !  albedo of ice under clear sky
+         palcnp                  !  albedo of ocean under clear sky
+    
+    !! * Local variables
+    INTEGER ::    &
+         ji, jj                   ! dummy loop indices
+    REAL(wp) ::   &
+         c1     = 0.05  , & ! constants values
+         c2     = 0.1   , &
+         albice = 0.50  , & !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
+         cgren  = 0.06  , & !  correction of the snow or ice albedo to take into account
+         !  effects of cloudiness (Grenfell & Perovich, 1984)
+         alphd  = 0.80  , & !  coefficients for linear interpolation used to compute
+         alphdi = 0.72  , & !  albedo between two extremes values (Pyane, 1972)
+         alphc  = 0.65  , &
+         zmue   = 0.4   , &     !  cosine of local solar altitude
+         zzero   = 0.0     ,  &
+         zone    = 1.0
+    
+    REAL(wp) ::   &
+         zmue14         , & !  zmue**1.4
+         zalbpsnm       , & !  albedo of ice under clear sky when snow is melting
+         zalbpsnf       , & !  albedo of ice under clear sky when snow is freezing
+         zalbpsn        , & !  albedo of snow/ice system when ice is coverd by snow
+         zalbpic        , & !  albedo of snow/ice system when ice is free of snow
+         zithsn         , & !  = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow)
+         zitmlsn        , & !  = 1 freezinz snow (sist >=rt0_snow) ; = 0 melting snow (sist<rt0_snow)
+         zihsc1         , & !  = 1 hsn <= c1 ; = 0 hsn > c1
+         zihsc2                   !  = 1 hsn >= c2 ; = 0 hsn < c2
+    REAL(wp), DIMENSION(jpi,jpj) ::  &
+         zalbfz         , & !  ( = alphdi for freezing ice ; = albice for melting ice )
+         zficeth                  !  function of ice thickness
+    LOGICAL , DIMENSION(jpi,jpj) ::  &
+         llmask
+    !!    to be included for without seaice
+    !!      REAL(wp), DIMENSION(jpi,jpj) ::   &  !:
+    !!      sist   ,   &  !: Sea-Ice Surface Temperature (Kelvin )
+    !!      hsnif  ,   &  !: Snow thickness
+    !!      hicif       !: Ice thickness
+    
+    !!---------------------------------------------------------------------
+    
+    !-------------------------
+    !  Computation of  zficeth
+    !--------------------------
+    
+    llmask = (hsnif == 0.0) .AND. ( sist >= rt0_ice )
+    WHERE ( llmask )   !  ice free of snow and melts
+       zalbfz = albice
+    ELSEWHERE
+       zalbfz = alphdi
+    END WHERE
+    
+    DO jj = 1, jpj
+       DO ji = 1, jpi
+          IF( hicif(ji,jj) > 1.5 ) THEN
+             zficeth(ji,jj) = zalbfz(ji,jj)
+          ELSEIF( hicif(ji,jj) > 1.0  .AND. hicif(ji,jj) <= 1.5 ) THEN
+             zficeth(ji,jj) = 0.472 + 2.0 * ( zalbfz(ji,jj) - 0.472 ) * ( hicif(ji,jj) - 1.0 )
+          ELSEIF( hicif(ji,jj) > 0.05 .AND. hicif(ji,jj) <= 1.0 ) THEN
+             zficeth(ji,jj) = 0.2467 + 0.7049 * hicif(ji,jj)                                &
+                  &                    - 0.8608 * hicif(ji,jj) * hicif(ji,jj)                 &
+                  &                    + 0.3812 * hicif(ji,jj) * hicif(ji,jj) * hicif (ji,jj)
+          ELSE
+             zficeth(ji,jj) = 0.1 + 3.6 * hicif(ji,jj)
+          ENDIF
+       END DO
+    END DO
+    
+    !-----------------------------------------------
+    !    Computation of the snow/ice albedo system
+    !-------------------------- ---------------------
+    
+    !    Albedo of snow-ice for clear sky.
+    !-----------------------------------------------
+     DO jj = 1, jpj
+        DO ji = 1, jpi
+           !  Case of ice covered by snow.
+           !  melting snow
+           zihsc1       = 1.0 - MAX ( zzero , SIGN ( zone , - ( hsnif(ji,jj) - c1 ) ) )
+           zalbpsnm     = ( 1.0 - zihsc1 ) & 
+                * ( zficeth(ji,jj) + hsnif(ji,jj) * ( alphd - zficeth(ji,jj) ) / c1 ) &
+                &                 + zihsc1   * alphd
+          !  freezing snow
+           zihsc2       = MAX ( zzero , SIGN ( zone , hsnif(ji,jj) - c2 ) )
+           zalbpsnf     = ( 1.0 - zihsc2 ) * & 
+                ( albice + hsnif(ji,jj) * ( alphc - albice ) / c2 ) &
+                &                 + zihsc2   * alphc
+           
+           zitmlsn      =  MAX ( zzero , SIGN ( zone , sist(ji,jj) - rt0_snow ) )
+           zalbpsn      =  zitmlsn * zalbpsnf + ( 1.0 - zitmlsn ) * zalbpsnm
+           
+           !  Case of ice free of snow.
+           zalbpic      = zficeth(ji,jj)
+           
+           ! albedo of the system
+           zithsn       = 1.0 - MAX ( zzero , SIGN ( zone , - hsnif(ji,jj) ) )
+           palbp(ji,jj) =  zithsn * zalbpsn + ( 1.0 - zithsn ) *  zalbpic
+        END DO
+     END DO
+    
+     !    Albedo of snow-ice for overcast sky.
+     !----------------------------------------------
+     palb(:,:)   = palbp(:,:) + cgren
+    
+     !--------------------------------------------
+     !    Computation of the albedo of the ocean
+     !-------------------------- -----------------
+    
+    
+     !  Parameterization of Briegled and Ramanathan, 1982
+     zmue14      = zmue**1.4
+     palcnp(:,:) = 0.05 / ( 1.1 * zmue14 + 0.15 )
+    
+     !  Parameterization of Kondratyev, 1969 and Payne, 1972
+     palcn(:,:)  = 0.06
+    
+   END SUBROUTINE flx_blk_albedo
+
+  
Index: /trunk/NEMO/OPA_SRC/cpl.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/cpl.F90	(revision 471)
+++ /trunk/NEMO/OPA_SRC/cpl.F90	(revision 472)
@@ -86,16 +86,7 @@
      
 #if defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_forced_daily
-      IF(lwp)WRITE(numout,cform_err)
-      IF(lwp)WRITE(numout,*) ' key_coupled and key_flx_bulk_... are incompatible'
-      nstop = nstop + 1
+      CALL ctl_stop( 'key_coupled and key_flx_bulk_... are incompatible' )
 #endif
  
-      IF(lwp)WRITE(numout,*)'     coupled simulation'
-      IF(lwp)WRITE(numout,*)'        unit ',numlhf,' receive atm fluxes'
-      IF(lwp)WRITE(numout,*)'        unit ',numlws,' receive windstress'
-      IF(lwp)WRITE(numout,*)'        unit ',numlts,' transfer sst'
-      IF(lwp)WRITE(numout,*)'        unit ',numlic,' transfer ice cover'
-
-
       CALL cpl_read           ! read the coupled mode information in namelist
 
Index: /trunk/NEMO/OPA_SRC/in_out_manager.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/in_out_manager.F90	(revision 471)
+++ /trunk/NEMO/OPA_SRC/in_out_manager.F90	(revision 472)
@@ -6,9 +6,4 @@
 
    PUBLIC
-   !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -32,25 +27,8 @@
       ninist = 0        ,    &  !: initial state output flag (0/1)
       nbench = 0                !: benchmark parameter (0/1)
-   !!----------------------------------------------------------------------
-   !!                          Run control  
-   !!----------------------------------------------------------------------
    
-   INTEGER ::                &  !:
-      nstop = 0 ,            &  !: e r r o r  flag (=number of reason for a
-      !                         !                   prematurely stop the run)
-      nwarn = 0                 !: w a r n i n g  flag (=number of warning
-      !                         !                       found during the run)
-
-   
-   CHARACTER (len=64) ::        &                                                    !:
-      cform_err="(/,' ===>>> : E R R O R',     /,'         ===========',/)"    ,   & !:
-      cform_war="(/,' ===>>> : W A R N I N G', /,'         ===============',/)"      !:
    !!----------------------------------------------------------------------
    !! output monitoring
    !! -----------------------------------
-
-   LOGICAL ::   &               !:
-      lwp                ,   &  !: boolean : true on the 1st processor only
-      lsp_area = .TRUE.         !: to make a control print over a specific area
 
    INTEGER ::                &  !:
@@ -80,27 +58,95 @@
       numisp     = 41 ,      &  !: logical unit for island statistics
       numgap     = 45 ,      &  !: logical unit for differences diagnostic
-      numwrs     = 46 ,      &  !: logical unit for output restart
-      numtdt     = 62 ,      &  !: logical unit for data temperature
-      numsdt     = 63 ,      &  !: logical unit for data salinity
-      numrnf     = 64 ,      &  !: logical unit for runoff data
-      numwso     = 71 ,      &  !: logical unit for 2d output write
-      numwvo     = 72 ,      &  !: logical unit for 3d output write
-      numsst     = 65 ,      &  !: logical unit for surface temperature data
-      numsss     = 66 ,      &  !: logical unit for surface salinity data
       numbol     = 67 ,      &  !: logical unit for "bol" diagnostics
       numptr     = 68 ,      &  !: logical unit for Poleward TRansports
-      numflo     = 69 ,      &  !: logical unit for drifting floats
+      numflo     = 69           !: logical unit for drifting floats
       !                         !: * coupled units
-      numlhf     = 71 ,      &  !: unit to transfer fluxes
-      numlws     = 72 ,      &  !: unit to transfer stress
-      numlts     = 73 ,      &  !: unit to transfer sst
-      numlic     = 74           !: unit to transfer ice cover
+
+   !!----------------------------------------------------------------------
+   !!                          Run control  
+   !!----------------------------------------------------------------------
+   
+   INTEGER ::                &  !:
+      nstop = 0 ,            &  !: e r r o r  flag (=number of reason for a
+      !                         !                   prematurely stop the run)
+      nwarn = 0                 !: w a r n i n g  flag (=number of warning
+      !                         !                       found during the run)
+
+   
+   CHARACTER(LEN=100) :: ctmp1, ctmp2, ctmp3    ! temporary character
+   CHARACTER (len=64) ::        &                                                    !:
+      cform_err="(/,' ===>>> : E R R O R',     /,'         ===========',/)"    ,   & !:
+      cform_war="(/,' ===>>> : W A R N I N G', /,'         ===============',/)"      !:
+   LOGICAL ::   &               !:
+      lwp                ,   &  !: boolean : true on the 1st processor only
+      lsp_area = .TRUE.         !: to make a control print over a specific area
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Header$ 
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !!----------------------------------------------------------------------
 
 
-   !! Contral/debugging
-   !! -----------------
-   REAL(wp) ::               &  !:
-      u_ctl, v_ctl,          &  !: sum of ua and va trend
-      t_ctl, s_ctl              !: sum of ta and sa trend
+CONTAINS
+
+
+   SUBROUTINE ctl_stop( cd1, cd2, cd3, cd4, cd5,   &
+      &                 cd6, cd7, cd8, cd9, cd10 )
+      !!-----------------------------------------------------------------------
+      !!                  ***  ROUTINE  stop_opa  ***
+      !!
+      !! ** Purpose : ???
+      !!
+      !!-----------------------------------------------------------------------
+      CHARACTER(len=*),INTENT(in),OPTIONAL ::  cd1, cd2, cd3, cd4, cd5, cd6, cd7, cd8, cd9, cd10
+      !!-----------------------------------------------------------------------
+      
+      nstop = nstop + 1 
+      IF(lwp) THEN
+         WRITE(numout,"(/,' ===>>> : E R R O R',     /,'         ===========',/)") 
+         IF( PRESENT(cd1 ) ) WRITE(numout,*) cd1
+         IF( PRESENT(cd2 ) ) WRITE(numout,*) cd2
+         IF( PRESENT(cd3 ) ) WRITE(numout,*) cd3
+         IF( PRESENT(cd4 ) ) WRITE(numout,*) cd4
+         IF( PRESENT(cd5 ) ) WRITE(numout,*) cd5
+         IF( PRESENT(cd6 ) ) WRITE(numout,*) cd6
+         IF( PRESENT(cd7 ) ) WRITE(numout,*) cd7
+         IF( PRESENT(cd8 ) ) WRITE(numout,*) cd8
+         IF( PRESENT(cd9 ) ) WRITE(numout,*) cd9
+         IF( PRESENT(cd10) ) WRITE(numout,*) cd10
+      ENDIF
+      CALL FLUSH(numout)
+
+   END SUBROUTINE ctl_stop
+
+
+   SUBROUTINE ctl_warn( cd1, cd2, cd3, cd4, cd5,   &
+      &                 cd6, cd7, cd8, cd9, cd10 )
+      !!-----------------------------------------------------------------------
+      !!                  ***  ROUTINE  stop_opa  ***
+      !!
+      !! ** Purpose : ???
+      !!
+      !!-----------------------------------------------------------------------
+      CHARACTER(len=*),INTENT(in),OPTIONAL ::  cd1, cd2, cd3, cd4, cd5, cd6, cd7, cd8, cd9, cd10
+      !!-----------------------------------------------------------------------
+      
+      nwarn = nwarn + 1 
+      IF(lwp) THEN
+         WRITE(numout,"(/,' ===>>> : W A R N I N G', /,'         ===============',/)") 
+         IF( PRESENT(cd1 ) ) WRITE(numout,*) cd1
+         IF( PRESENT(cd2 ) ) WRITE(numout,*) cd2
+         IF( PRESENT(cd3 ) ) WRITE(numout,*) cd3
+         IF( PRESENT(cd4 ) ) WRITE(numout,*) cd4
+         IF( PRESENT(cd5 ) ) WRITE(numout,*) cd5
+         IF( PRESENT(cd6 ) ) WRITE(numout,*) cd6
+         IF( PRESENT(cd7 ) ) WRITE(numout,*) cd7
+         IF( PRESENT(cd8 ) ) WRITE(numout,*) cd8
+         IF( PRESENT(cd9 ) ) WRITE(numout,*) cd9
+         IF( PRESENT(cd10) ) WRITE(numout,*) cd10
+      ENDIF
+      CALL FLUSH(numout)
+
+   END SUBROUTINE ctl_warn
 
 END MODULE in_out_manager
Index: /trunk/NEMO/OPA_SRC/iom.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/iom.F90	(revision 472)
+++ /trunk/NEMO/OPA_SRC/iom.F90	(revision 472)
@@ -0,0 +1,733 @@
+MODULE iom
+   !!=====================================================================
+   !!                    ***  MODULE  iom ***
+   !!
+   !! Input/Output manager :  Library to read input files
+   !!
+   !! Ongoing work : This code is here to help discussions about I/O
+   !!                library in the NEMO system
+   !!====================================================================
+   !!--------------------------------------------------------------------
+   !!   iom_open       : open a file read only
+   !!   iom_close      : close a file or all files opened by iom
+   !!   iom_get        : read a field : interface to several routines
+   !!   iom_varid      : get the id of a variable in a file
+   !!   iom_get_gblatt : ???
+   !!--------------------------------------------------------------------
+   !! History :  9.0  ! 05 12  (J. Belier) Original code
+   !!            9.0  ! 06 02  (S. Masson) Adaptation to NEMO
+   !!--------------------------------------------------------------------
+   !! * Modules used
+   USE in_out_manager  ! I/O manager
+   USE dom_oce         ! ocean space and time domain
+   USE lbclnk          ! ???
+   USE ioipsl          ! ???
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC :: iom_open, iom_close, iom_get, iom_varid, iom_get_gblatt
+
+   !! * Interfaces
+   INTERFACE iom_get
+      MODULE PROCEDURE iom_get_r_1d, iom_get_r_2d, iom_get_r_3d
+   END INTERFACE
+
+   !! * Share module variables
+   INTEGER, PARAMETER, PUBLIC ::        &  !:
+      jpdom_data                  = 1,  &  !: ( 1  :jpidta, 1  :jpjdta)
+      jpdom_global                = 2,  &  !: ( 1  :jpiglo, 1  :jpjglo)
+      jpdom_local                 = 3,  &  !: One of the 3 following cases
+      jpdom_local_full            = 4,  &  !: ( 1  :jpi   , 1  :jpi   )
+      jpdom_local_noextra         = 5,  &  !: ( 1  :nlci  , 1  :nlcj  )
+      jpdom_local_noovlap         = 6,  &  !: (nldi:nlei  ,nldj:nlej  )
+      jpdom_unknown               = 7      !: No dimension checking
+
+   !! * Module variables
+   INTEGER, PARAMETER ::    &
+      jpmax_vars    = 50,   &  ! maximum number of variables in one file
+      jpmax_dims    =  5,   &  ! maximum number of dimensions for one variable
+      jpmax_digits  =  5       ! maximum number of digits in the file name to reference the cpu number
+ 
+   INTEGER :: iom_init = 0
+
+   TYPE :: flio_file
+      CHARACTER(LEN=240)                        ::   name     ! name of the file
+      INTEGER                                   ::   iopen    ! 1/0 is the file is open/not open
+      INTEGER                                   ::   nvars    ! number of identified varibles in the file
+      INTEGER                                   ::   iduld    ! id of the unlimited dimension
+      CHARACTER(LEN=16), DIMENSION(jpmax_vars)  ::   cn_var   ! names of the variables
+      INTEGER, DIMENSION(jpmax_vars)            ::   ndims    ! number of dimensions of the variables
+      LOGICAL, DIMENSION(jpmax_vars)            ::   luld     ! variable including unlimited dimension
+      INTEGER, DIMENSION(jpmax_dims,jpmax_vars) ::   dimsz    ! size of the dimensions of the variables
+      REAL(kind=wp), DIMENSION(jpmax_vars)      ::   scf      ! scale_factor of the variables
+      REAL(kind=wp), DIMENSION(jpmax_vars)      ::   ofs      ! add_offset of the variables
+   END TYPE flio_file
+   TYPE(flio_file), DIMENSION(flio_max_files)   :: iom_file ! array containing the info for all opened files
+
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2006)
+   !! $Header$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE iom_open( cdname, knumfl, ldimg )
+      !!---------------------------------------------------------------------
+      !!                   ***  SUBROUTINE  iom_open  ***
+      !!
+      !! ** Purpose :  open an input file read only (return 0 if not found)
+      !!
+      !! ** Method :
+      !!
+      !!---------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in )  ::   cdname   ! File name
+      INTEGER, INTENT(out)           ::   knumfl   ! Identifier of the opened file
+      LOGICAL, INTENT(in ), OPTIONAL ::   ldimg    ! flg to specify that we use dimg format
+
+      CHARACTER(LEN=100) :: clname   ! the name of the file based on cdname [[+clcpu]+clcpu]
+      CHARACTER(LEN=10) :: clsuffix  ! ".nc" or ".dimg"
+      CHARACTER(LEN=10) :: clcpu     ! the cpu number (max jpmax_digits digits)
+      LOGICAL :: llok                ! check the existence 
+      INTEGER :: icnt                ! counter for digits in clcpu (max = jpmax_digits)
+      !---------------------------------------------------------------------
+
+      ! find the file
+      ! =============
+      clname   = trim(cdname)
+#if defined key_agrif
+      if ( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
+#endif                
+      clsuffix = ".nc"
+      IF( PRESENT(ldimg) ) THEN
+         IF ( ldimg ) clsuffix = ".dimg"
+      ENDIF
+      !
+      INQUIRE( FILE = clname, EXIST = llok )
+      IF( .NOT.llok ) THEN                         ! try to complete the name with the suffix only
+         clname = TRIM(cdname)//TRIM(clsuffix)
+         INQUIRE( FILE = clname, EXIST = llok )
+         IF( .NOT.llok ) THEN                      ! try to complete the name with both cpu number and suffix
+            WRITE(clcpu,*) narea-1
+            clcpu  = trim(adjustl(clcpu))
+            clname = trim(cdname)//"_"
+            icnt = 0
+            INQUIRE( FILE = trim(clname)//trim(clcpu)//trim(clsuffix), EXIST = llok ) 
+            DO WHILE( .NOT.llok .AND. icnt < jpmax_digits )   ! we try fifferent formats for the cpu number by adding 0
+               clname = trim(clname)//"0"
+               INQUIRE( FILE = trim(clname)//trim(clcpu)//trim(clsuffix), EXIST = llok )
+               icnt = icnt + 1
+            END DO
+            IF( .NOT.llok ) THEN                   ! no way to find the files...
+               CALL ctl_stop( 'iom_open: file '//trim(clname)//'... not found' )
+            ENDIF
+            clname = trim(clname)//trim(clcpu)//trim(clsuffix)
+         ENDIF
+      ENDIF
+
+      ! Open the file
+      ! =============
+      IF( llok ) THEN                         
+         IF (lwp) WRITE(numout,*) 'iom_open ~~~  open file: '//trim(clname)
+         CALL flioopfd( trim(clname), knumfl )
+         IF( iom_init == 0 ) THEN
+            iom_file(:)%iopen = 0
+            iom_init = 1
+         ENDIF
+         iom_file(knumfl)%iopen      = 1
+         iom_file(knumfl)%name       = TRIM(clname)
+         iom_file(knumfl)%nvars      = 0
+         iom_file(knumfl)%ndims(:)   = 0
+         iom_file(knumfl)%luld(:)    = .FALSE.
+         iom_file(knumfl)%dimsz(:,:) = 0
+         ! does the file contain time axis (that must be unlimitted) ?
+         CALL flioinqf( knumfl, id_uld = iom_file(knumfl)%iduld )
+      ELSE
+         knumfl = 0
+      ENDIF
+  
+   END SUBROUTINE iom_open
+
+
+   SUBROUTINE iom_close( knumfl )
+      !!--------------------------------------------------------------------
+      !!                   ***  SUBROUTINE  iom_close  ***
+      !!
+      !! ** Purpose : close an input file, or all files opened by iom
+      !!
+      !! ** Method :
+      !!
+      !!--------------------------------------------------------------------
+      INTEGER,INTENT(in), OPTIONAL ::   knumfl   ! Identifier of the file to be closed
+      !                                          ! If this argument is not present,
+      !                                          ! all the files opened by iom are closed.
+
+      INTEGER           ::   jf            ! dummy loop indices
+      INTEGER           ::   i_s, i_e      ! temporary integer
+      !---------------------------------------------------------------------
+
+      IF( PRESENT(knumfl) ) THEN
+         i_s = knumfl
+         i_e = knumfl
+      ELSE
+         i_s = 1
+         i_e = flio_max_files
+      ENDIF
+      IF ( i_s > 0 ) THEN
+         DO jf = i_s, i_e
+            IF( iom_file(jf)%iopen > 0 ) THEN
+               CALL flioclo( jf )
+               iom_file(jf)%iopen      = 0
+               iom_file(jf)%name       = 'NONE'
+               iom_file(jf)%nvars      = 0
+               iom_file(jf)%iduld      = 0
+               iom_file(jf)%ndims(:)   = 0
+               iom_file(jf)%luld(:)    = .FALSE.
+               iom_file(jf)%dimsz(:,:) = 0
+            ELSEIF( PRESENT(knumfl) ) THEN
+               WRITE(ctmp1,*) '--->',  knumfl
+               CALL ctl_stop( 'iom_close: Invalid file identifier', ctmp1 )
+            ENDIF
+         END DO
+      ENDIF
+          
+   END SUBROUTINE iom_close
+ 
+
+   !!----------------------------------------------------------------------
+   !!                   INTERFACE iom_u_getv
+   !!----------------------------------------------------------------------
+   SUBROUTINE iom_get_r_1d( knumfl, kdom , cdvar , pvar  ,   &
+      &                             ktime, kstart, kcount )
+      INTEGER               , INTENT(in )           ::   knumfl    ! Identifier of the file
+      INTEGER               , INTENT(in )           ::   kdom      ! Type of domain to be read
+      CHARACTER(len=*)      , INTENT(in )           ::   cdvar     ! Name of the variable
+      REAL(wp), DIMENSION(:), INTENT(out)           ::   pvar      ! read field
+      INTEGER               , INTENT(in ) ,OPTIONAL ::   ktime     ! record number
+      INTEGER , DIMENSION(:), INTENT(in ), OPTIONAL ::   kstart    ! start position of the reading in each axis 
+      INTEGER , DIMENSION(:), INTENT(in ), OPTIONAL ::   kcount    ! number of points to be read in each axis
+!
+      CHARACTER(LEN=100) :: clinfo                    ! info character
+!
+      clinfo = 'iom_get_r_1d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
+      IF( PRESENT(kstart) ) THEN
+         IF ( SIZE(kstart) /= 1 ) CALL ctl_stop( trim(clinfo), 'kstart must be a 1 element vector' )
+      ENDIF
+      IF( PRESENT(kcount) ) THEN
+         IF ( SIZE(kcount) /= 1 ) CALL ctl_stop( trim(clinfo), 'kcount must be a 1 element vector' )
+      ENDIF
+      IF ( knumfl > 0 ) CALL iom_u_getv( knumfl, kdom       , cdvar        , pv_r1d=pvar,   &
+           &                                     ktime=ktime, kstart=kstart, kcount=kcount )
+   END SUBROUTINE iom_get_r_1d
+   SUBROUTINE iom_get_r_2d( knumfl, kdom , cdvar , pvar  ,   &
+      &                             ktime, kstart, kcount )
+      INTEGER,INTENT(in) :: knumfl
+      INTEGER,INTENT(in) :: kdom
+      CHARACTER(len=*),INTENT(in) :: cdvar
+      REAL(wp),INTENT(out),DIMENSION(:,:) :: pvar
+      INTEGER,INTENT(in),OPTIONAL :: ktime
+      INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kstart
+      INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kcount
+!
+      CHARACTER(LEN=100) :: clinfo                    ! info character
+!
+      clinfo = 'iom_get_r_2d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
+      IF( PRESENT(kstart) ) THEN
+         IF ( size(kstart) /= 2 ) CALL ctl_stop(trim(clinfo), 'kstart must be a 2 element vector')
+      ENDIF
+      IF( PRESENT(kcount) ) THEN
+         IF ( size(kcount) /= 2 ) CALL ctl_stop(trim(clinfo), 'kcount must be a 2 element vector')
+      ENDIF
+      IF ( knumfl > 0 ) CALL iom_u_getv( knumfl, kdom       , cdvar        , pv_r2d=pvar,   &
+         &                                       ktime=ktime, kstart=kstart, kcount=kcount )
+   END SUBROUTINE iom_get_r_2d
+   SUBROUTINE iom_get_r_3d( knumfl, kdom , cdvar , pvar  ,   &
+      &                             ktime, kstart, kcount )
+      INTEGER,INTENT(in) :: knumfl
+      INTEGER,INTENT(in) :: kdom
+      CHARACTER(len=*),INTENT(in) :: cdvar
+      REAL(wp),INTENT(out),DIMENSION(:,:,:) :: pvar
+      INTEGER,INTENT(in),OPTIONAL :: ktime
+      INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kstart
+      INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kcount
+!
+      CHARACTER(LEN=100) :: clinfo                    ! info character
+!
+      clinfo = 'iom_get_r_3d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
+      IF ( PRESENT(kstart) ) THEN
+         IF ( size(kstart) /= 3 ) CALL ctl_stop(trim(clinfo), 'kstart must be a 3 element vector')
+      ENDIF
+      IF ( PRESENT(kcount) ) THEN
+         IF ( size(kcount) /= 3 ) CALL ctl_stop(trim(clinfo), 'kcount must be a 3 element vector')
+      ENDIF
+      IF ( knumfl > 0 ) CALL iom_u_getv( knumfl, kdom       , cdvar        , pv_r3d=pvar,   &
+	     &                                   ktime=ktime, kstart=kstart, kcount=kcount )
+   END SUBROUTINE iom_get_r_3d
+   !!----------------------------------------------------------------------
+
+
+   SUBROUTINE iom_u_getv( knumfl, kdom  , cdvar , &
+        &                   pv_r1d, pv_r2d, pv_r3d, &
+        &                   ktime , kstart, kcount )
+     !!-----------------------------------------------------------------------
+     !!                  ***  ROUTINE  iom_u_getv  ***
+     !!
+     !! ** Purpose : read a 1D/2D/3D variable
+     !!
+     !! ** Method : read ONE time step at each CALL
+     !!
+     !!-----------------------------------------------------------------------
+     INTEGER,                      INTENT(in )           ::   knumfl     ! Identifier of the file
+     INTEGER,                      INTENT(in )           ::   kdom       ! Type of domain to be read
+     CHARACTER(len=*),             INTENT(in )           ::   cdvar      ! Name of the variable
+     REAL(wp), DIMENSION(:)      , INTENT(out), OPTIONAL ::   pv_r1d     ! read field (1D case)
+     REAL(wp), DIMENSION(:,:)    , INTENT(out), OPTIONAL ::   pv_r2d     ! read field (2D case)
+     REAL(wp), DIMENSION(:,:,:)  , INTENT(out), OPTIONAL ::   pv_r3d     ! read field (3D case)
+     INTEGER                     , INTENT(in ), OPTIONAL ::   ktime      ! record number
+     INTEGER , DIMENSION(:)      , INTENT(in ), OPTIONAL ::   kstart     ! start position of the reading in each axis 
+     INTEGER , DIMENSION(:)      , INTENT(in ), OPTIONAL ::   kcount     ! number of points to be read in each axis
+
+     INTEGER                        :: jl            ! loop on number of dimension 
+     INTEGER                        :: idom,    &    ! type of domain
+          &                            idvar,   &    ! id of the variable
+          &                            inbdim,  &    ! number of dimensions of the variable
+          &                            idmspc,  &    ! number of spatial dimensions 
+          &                            itime,   &    ! record number
+          &                            istop         ! temporary value of nstop
+     INTEGER, DIMENSION(jpmax_dims) :: istart,  &    ! starting point to read for each axis
+          &                            icnt,    &    ! number of value to read along each axis 
+          &                            idimsz        ! size of the dimensions of the variable
+     REAL(wp)                       :: zscf, zofs    ! sacle_factor and add_offset
+     INTEGER                        ::  itmp         ! temporary integer
+     CHARACTER(LEN=100)             :: clinfo        ! info character
+     !---------------------------------------------------------------------
+     clinfo = 'iom_u_getv, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
+     ! local definition of the domain ?
+     idom = kdom
+     ! check kcount and kstart optionals parameters...
+     IF( PRESENT(kcount) .AND. (.NOT. PRESENT(kstart)) ) &
+          CALL ctl_stop( trim(clinfo), 'kcount present needs kstart present' )
+     IF( PRESENT(kstart) .AND. (.NOT. PRESENT(kcount)) ) &
+          CALL ctl_stop( trim(clinfo), 'kstart present needs kcount present' )
+     IF( PRESENT(kstart) .AND. idom /= jpdom_unknown ) &
+          CALL ctl_stop( trim(clinfo), 'kstart present needs kdom = jpdom_unknown' )
+
+     ! Search for the variable in the data base (eventually actualize data)
+     !-
+     istop = nstop
+     idvar = iom_varid( knumfl, cdvar )
+     !
+     IF ( idvar > 0 ) THEN
+        ! to write iom_file(knumfl)%dimsz in a shorter way !
+        idimsz(:) = iom_file(knumfl)%dimsz(:, idvar) 
+        inbdim = iom_file(knumfl)%ndims(idvar)! number of dimensions in the file
+        idmspc = inbdim ! number of spatial dimensions in the file
+        IF( iom_file(knumfl)%luld(idvar) ) idmspc = inbdim - 1
+        IF( idmspc > 3 ) CALL ctl_stop(trim(clinfo), 'the file has more than 3 spatial dimensions',   &
+             &                    'this case is not coded...') 
+        ! Identify the domain in case of jpdom_local definition
+        !-
+        IF( idom == jpdom_local ) THEN
+           IF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN 
+              idom = jpdom_local_full
+           ELSEIF( idimsz(1) == nlci .AND. idimsz(2) == nlcj ) THEN 
+              idom = jpdom_local_noextra
+           ELSEIF( idimsz(1) == (nlei - nldi + 1) .AND. idimsz(2) == (nlej - nldj + 1) ) THEN 
+              idom = jpdom_local_noovlap
+           ELSE 
+              CALL ctl_stop( trim(clinfo), 'impossible to identify the local domain' )
+           ENDIF
+        ENDIF
+
+        ! definition of istart and icnt
+        !-
+        ! initializations
+        istart(:) = 1
+        icnt(:) = 1
+        itime = 1
+        IF( PRESENT(ktime) ) itime = ktime
+        !
+        IF( PRESENT(pv_r1d) ) THEN
+           IF( idmspc == 1 ) THEN 
+              ! data is 1d array (+ maybe a temporal dimension)
+              IF( PRESENT(kstart) ) THEN
+                 istart(1:2) = (/ kstart(1), itime /)
+                 icnt(1) = kcount(1)
+              ELSE
+                 IF( kdom == jpdom_unknown ) THEN
+                    istart(2) = itime
+                    icnt(1) = idimsz(1)
+                 ELSE 
+                    CALL ctl_stop( trim(clinfo), 'case not coded...You must use jpdom_unknown' )
+                 ENDIF
+              ENDIF
+           ELSE
+              CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 1D array,', &
+                   &         'we do not accept data with more than 1 spatial dimension',     &
+                   &         'Use ncwa -a to suppress the unnecessary dimensions')
+           ENDIF
+        ELSEIF( PRESENT(pv_r2d) ) THEN
+           SELECT CASE (idmspc)
+           CASE (1)
+              CALL ctl_stop( trim(clinfo), 'the file has only 1 spatial dimension',  &
+                   &         'it is impossible to read a 2d array from this file...')
+           CASE (2)
+              ! data is 2d array (+ maybe a temporal dimension)
+              IF ( PRESENT(kstart) ) THEN
+                 istart(1:3) = (/ kstart(1:2), itime /)
+                 icnt(1:2) = kcount(1:2)
+              ELSE
+                 IF( kdom == jpdom_unknown ) THEN
+                    istart(3) = itime
+                    icnt(1:2) = idimsz(1:2)
+                 ELSE 
+                    IF( idom == jpdom_data ) THEN
+                       istart(1:3) = (/ mig(1), mjg(1), itime /)
+                    ELSEIF( idom == jpdom_global ) THEN
+                       istart(1:3) = (/ nimpp, njmpp, itime /)
+                    ENDIF
+                    ! we do not read the overlap -> we start to read at nldi, nldj
+                    IF( idom /= jpdom_local_noovlap )   istart(1:2) = istart(1:2) + (/ nldi - 1, nldj - 1 /)
+                    ! we do not read the overlap and the extra-halos  
+                    ! -> from nldi to nlei and from nldj to nlej 
+                    icnt(1:2) = (/ nlei - nldi + 1, nlej - nldj + 1 /)
+                 ENDIF
+              ENDIF
+           CASE DEFAULT
+              IF ( itime == 1 .AND. idimsz(3) == 1 .AND. idmspc == 3 ) THEN
+                 CALL ctl_warn( trim(clinfo), '2D array but 3 spatial dimensions for the data...', &
+                      &         'As the size of the z dimension is 1 and as we try to read the first reacord, ',     &
+                      &         'we accept this case even if there is a possible mix-up between z and time dimension...')           
+                 IF ( PRESENT(kstart) ) THEN
+                    istart(1:2) = kstart(1:2)
+                    icnt(1:2) = kcount(1:2)
+                 ELSE
+                    IF( kdom == jpdom_unknown ) THEN
+                       icnt(1:2) = idimsz(1:2)
+                    ELSE 
+                       IF( idom == jpdom_data ) THEN
+                          istart(1:2) = (/ mig(1), mjg(1) /)
+                       ELSEIF( idom == jpdom_global ) THEN
+                          istart(1:2) = (/ nimpp, njmpp /)
+                       ENDIF
+                       ! we do not read the overlap -> we start to read at nldi, nldj
+                       IF( idom /= jpdom_local_noovlap )   istart(1:2) = istart(1:2) + (/ nldi - 1, nldj - 1 /)
+                       ! we do not read the overlap and the extra-halos  
+                       ! -> from nldi to nlei and from nldj to nlej 
+                       icnt(1:2) = (/ nlei - nldi + 1, nlej - nldj + 1 /)
+                    ENDIF
+                 ENDIF
+              ELSE
+                 CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 2D array,', &
+                      &         'we do not accept data with more than 2 spatial dimension',     &
+                      &         'Use ncwa -a to suppress the unnecessary dimensions')           
+              ENDIF
+           END SELECT
+        ELSEIF( PRESENT(pv_r3d) ) THEN
+           SELECT CASE (idmspc)
+           CASE (1)
+              CALL ctl_stop(trim(clinfo), 'the file has only 1 spatial dimension',   &
+                   &        'it is impossible to read a 3d array from this file...')
+           CASE (2)
+              CALL ctl_stop(trim(clinfo), 'the file has only 2 spatial dimension',   &
+                   &        'it is impossible to read a 3d array from this file...')
+           CASE (3)
+              ! data is 3d array (+ maybe a temporal dimension)
+              IF( PRESENT(kstart) ) THEN
+                 istart(1:4) = (/ kstart(1:3), itime /)
+                 icnt(1:3) = kcount(1:3)
+              ELSE
+                 IF( kdom == jpdom_unknown ) THEN
+                    istart(4) = itime
+                    icnt(1:3) = idimsz(1:3)
+                 ELSE 
+                    IF( idom == jpdom_data ) THEN
+                       istart(1:4) = (/ mig(1), mjg(1), 1, itime /)
+                    ELSEIF( idom == jpdom_global ) THEN
+                       istart(1:4) = (/ nimpp, njmpp, 1, itime /)
+                    ENDIF
+                    ! we do not read the overlap -> we start to read at nldi, nldj
+                    IF( idom /= jpdom_local_noovlap )   istart(1:2) = istart(1:2) + (/ nldi - 1, nldj - 1 /)
+                    ! we do not read the overlap and the extra-halos
+                    ! -> from nldi to nlei and from nldj to nlej 
+                    icnt(1:2) = (/ nlei - nldi + 1, nlej - nldj + 1 /)
+                    IF( idom == jpdom_data ) THEN 
+                       icnt(3) = jpkdta
+                    ELSE
+                       icnt(3) = jpk
+                    ENDIF
+                 ENDIF
+              ENDIF
+           CASE DEFAULT
+              CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 3D array,', &
+                   &         'we do not accept data with more than 3 spatial dimension',     &
+                   &         'Use ncwa -a to suppress the unnecessary dimensions')           
+           END SELECT
+        ENDIF
+
+        ! check that istart and icnt can be used with this file
+        !-
+        DO jl = 1, jpmax_dims
+           itmp = istart(jl)+icnt(jl)-1
+           IF( (itmp) > idimsz(jl) .AND. idimsz(jl) /= 0 ) THEN
+              WRITE(ctmp1,*) '(istart(',jl,') + icnt(',jl,') - 1) = ', itmp
+              WRITE(ctmp2,*) ' is larger than idimsz(',jl,'): ', idimsz(jl)
+              CALL ctl_stop( trim(clinfo), 'start and count too big regarding to the size of the data, ', ctmp1, ctmp2 )     
+           ENDIF
+        END DO
+
+        ! check that istart and icnt match the input array
+        !-     
+        IF( PRESENT(pv_r1d) ) THEN
+           itmp = size(pv_r1d)
+           WRITE(ctmp1,*) 'size(pv_r1d): ', itmp, ' /= icnt(1): ', icnt(1)
+           IF( itmp /= icnt(1) ) CALL ctl_stop( trim(clinfo), ctmp1 )
+        ELSEIF( PRESENT(pv_r2d) ) THEN
+           DO jl = 1, 2
+              IF( idom == jpdom_unknown ) THEN
+                 itmp = size(pv_r2d, jl)
+                 WRITE(ctmp1,*) 'size(pv_r2d, ',jl,'): ',itmp,' /= icnt(',jl,'): ',icnt(jl)
+              ELSE 
+                 itmp = size(pv_r2d(nldi:nlei,nldj:nlej), jl)
+                 WRITE(ctmp1,*) 'size(pv_r2d(nldi:nlei,nldj:nlej), ',jl,'): ',itmp,' /= icnt(',jl,'): ',icnt(jl)
+              ENDIF
+              IF( itmp /= icnt(jl) ) CALL ctl_stop( trim(clinfo), ctmp1 )
+           END DO
+        ELSEIF( PRESENT(pv_r3d) ) THEN
+           DO jl = 1, 3
+              IF( idom == jpdom_unknown ) THEN
+                 itmp = size(pv_r3d, jl)
+                 WRITE(ctmp1,*) 'size(pv_r3d, ',jl,'): ',itmp,' /= icnt(',jl,'):', icnt(jl)
+              ELSE 
+                 itmp = size(pv_r3d(nldi:nlei,nldj:nlej,:), jl)
+                 WRITE(ctmp1,*) 'size(pv_r3d(nldi:nlei,nldj:nlej), ',jl,'): ',itmp,' /= icnt(',jl,'): ',icnt(jl)
+              ENDIF
+              IF( itmp /= icnt(jl) ) CALL ctl_stop( trim(clinfo), ctmp1 )
+           END DO
+        ENDIF
+     ENDIF
+
+     ! read the data
+     !-     
+     IF( istop == nstop) THEN ! no additional errors until this point...
+        !
+        istop = nstop
+        !
+        zscf = iom_file(knumfl)%scf(idvar)      ! scale factor
+        zofs = iom_file(knumfl)%ofs(idvar)      ! offset
+        !
+        IF( PRESENT(pv_r1d) ) THEN
+           CALL fliogetv( knumfl, cdvar, pv_r1d(:), start=istart(1:inbdim), count=icnt(1:inbdim) )
+           !--- Apply scale_factor and offset
+           IF( zscf /= 1. ) pv_r1d(:) = pv_r1d(:) * zscf 
+           IF( zofs /= 0. ) pv_r1d(:) = pv_r1d(:) + zofs
+        ELSEIF( PRESENT(pv_r2d) ) THEN
+           IF( idom /= jpdom_unknown ) THEN
+              CALL fliogetv( knumfl, cdvar, pv_r2d(nldi:nlei,nldj:nlej), start=istart(1:inbdim), count=icnt(1:inbdim) )
+              !--- Apply scale_factor and offset
+              IF (zscf /= 1.) pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) * zscf
+              IF (zofs /= 0.) pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) + zofs
+              !--- Fill the overlap areas and extra hallows (mpp)
+              CALL lbc_lnk (pv_r2d,'Z',-999.,'no0')
+           ELSE
+              CALL fliogetv( knumfl, cdvar, pv_r2d(:,:), start=istart(1:inbdim), count=icnt(1:inbdim) )
+              !--- Apply scale_factor and offset
+              IF (zscf /= 1.) pv_r2d(:,:) = pv_r2d(:,:) * zscf
+              IF (zofs /= 0.) pv_r2d(:,:) = pv_r2d(:,:) + zofs
+           ENDIF
+        ELSEIF( PRESENT(pv_r3d) ) THEN
+           IF( idom /= jpdom_unknown ) THEN
+              CALL fliogetv( knumfl, cdvar, pv_r3d(nldi:nlei,nldj:nlej,:), start=istart(1:inbdim), count=icnt(1:inbdim) )
+              !--- Apply scale_factor and offset
+              IF( zscf /= 1. ) pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) * zscf
+              IF( zofs /= 0. ) pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) + zofs
+              !--- Fill the overlap areas and extra hallows (mpp)
+              IF( jpk == jpkdta )   CALL lbc_lnk( pv_r3d,'Z',-999.,'no0' ) ! this if could be removed with the new lbc_lnk ...
+           ELSE
+              CALL fliogetv( knumfl, cdvar, pv_r3d(:,:,:), start=istart(1:inbdim), count=icnt(1:inbdim) )
+              !--- Apply scale_factor and offset
+              IF (zscf /= 1.)   pv_r3d(:,:,:) = pv_r3d(:,:,:) * zscf
+              IF (zofs /= 0.)   pv_r3d(:,:,:) = pv_r3d(:,:,:) + zofs
+           ENDIF
+        ENDIF
+        !
+        IF( istop == nstop .AND. lwp )  &
+             &  WRITE(numout,*) ' read '//trim(cdvar)//' in '//trim(iom_file(knumfl)%name)//' ok'
+     ENDIF
+     !
+   END SUBROUTINE iom_u_getv
+
+   
+   SUBROUTINE iom_gettime( knumfl, cdvar, ptime )
+     !!--------------------------------------------------------------------
+     !!                   ***  SUBROUTINE  iom_close  ***
+     !!
+     !! ** Purpose : read the time axis cdvar in the file 
+     !!
+     !! ** Method :
+     !!
+     !!--------------------------------------------------------------------
+     INTEGER               , INTENT(in)  ::   knumfl   ! Identifier of the file to be closed
+     CHARACTER(len=*)      , INTENT(in)  ::   cdvar    ! time axis name
+     REAL(wp), DIMENSION(:), INTENT(out) ::   ptime    ! the time axis
+
+     INTEGER           :: idvar    ! id of the variable
+     CHARACTER(LEN=100) :: clinfo                    ! info character
+     !---------------------------------------------------------------------
+     clinfo = 'iom_gettime, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
+     idvar = iom_varid( knumfl, cdvar )
+     !
+     ptime(:) = 0. ! default definition
+     IF ( idvar > 0 ) THEN
+        IF ( iom_file(knumfl)%ndims(idvar) == 1 ) THEN
+           IF ( iom_file(knumfl)%luld(idvar) ) THEN
+              IF ( iom_file(knumfl)%dimsz(1,idvar) == size(ptime) ) THEN
+                 CALL fliogetv( knumfl, cdvar, ptime(:), start=(/ 1 /), &
+                      &                                   count=(/ iom_file(knumfl)%dimsz(1,idvar) /) )
+              ELSE
+                 WRITE(ctmp1,*) 'error with the size of ptime ',size(ptime),iom_file(knumfl)%dimsz(1,idvar)
+                 CALL ctl_stop( trim(clinfo), trim(ctmp1) )
+              ENDIF
+           ELSE
+              CALL ctl_stop( trim(clinfo), 'variable dimension is not unlimited... use iom_get' )
+           ENDIF
+        ELSE
+           CALL ctl_stop( trim(clinfo), 'the variable has more than 1 dimension' )
+        ENDIF
+     ELSE
+        CALL ctl_stop( trim(clinfo), 'variable not found in '//iom_file(knumfl)%name )
+     ENDIF
+
+   END SUBROUTINE iom_gettime
+    
+    
+    
+   FUNCTION iom_varid ( knumfl, cdvar, kdimsz )  
+     !!-----------------------------------------------------------------------
+     !!                  ***  FUNCTION  iom_varid  ***
+     !!
+     !! ** Purpose : get the id of a variable in a file (return 0 if not found)
+     !!
+     !! ** Method : ???
+     !!
+     !!-----------------------------------------------------------------------
+     INTEGER              , INTENT(in)            ::   knumfl   ! file Identifier
+     CHARACTER(len=*)     , INTENT(in)            ::   cdvar    ! name of the variable
+     INTEGER, DIMENSION(:), INTENT(out), OPTIONAL ::   kdimsz   ! size of the dimensions
+     !
+     INTEGER                        :: iom_varid, idvar, i_nvd, ji
+     INTEGER, DIMENSION(jpmax_dims) :: idimid
+     LOGICAL                        :: ll_fnd
+     CHARACTER(LEN=100)             :: clinfo                   ! info character
+     !!-----------------------------------------------------------------------
+     clinfo = 'iom_varid, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
+     iom_varid = 0                         ! default definition
+     IF ( PRESENT(kdimsz) ) kdimsz(:) = 0  ! default definition
+     !
+     IF ( knumfl > 0 ) THEN
+        IF( iom_file(knumfl)%iopen == 0 ) THEN 
+           CALL ctl_stop( trim(clinfo), 'the file is not open' )
+        ELSE
+           !
+           ll_fnd  = .FALSE.
+           idvar = 0
+           !
+           DO WHILE ( .NOT.ll_fnd .AND. idvar < iom_file(knumfl)%nvars )
+              idvar = idvar + 1
+              ll_fnd  = ( TRIM(cdvar) == TRIM(iom_file(knumfl)%cn_var(idvar)) )
+           END DO
+           !
+           IF( .NOT.ll_fnd ) THEN
+              idvar = idvar + 1
+              IF( idvar <= jpmax_vars ) THEN
+                 CALL flioinqv( knumfl, cdvar, ll_fnd, nb_dims = i_nvd )
+                 IF( ll_fnd ) THEN
+                    IF( i_nvd <= jpmax_dims ) THEN
+                       iom_file(knumfl)%nvars           = idvar
+                       iom_file(knumfl)%cn_var(idvar) = trim(cdvar)
+                       iom_file(knumfl)%ndims(idvar)  = i_nvd
+                       CALL flioinqv( knumfl, cdvar, ll_fnd, &
+                            &         len_dims = iom_file(knumfl)%dimsz(1:i_nvd,idvar), &
+                            &         id_dims  = idimid(1:i_nvd) )
+                       DO ji = 1, i_nvd
+                          IF ( idimid(ji) == iom_file(knumfl)%iduld ) iom_file(knumfl)%luld(idvar) = .TRUE.
+                       END DO
+                       !----------
+                       !---------- Deal with scale_factor and offset
+                       CALL flioinqa( knumfl, cdvar, 'scale_factor', ll_fnd )
+                       IF (ll_fnd) THEN
+                          CALL fliogeta( knumfl, cdvar, 'scale_factor', iom_file(knumfl)%scf(idvar) )
+                       ELSE
+                          iom_file(knumfl)%scf(idvar) = 1.
+                       END IF
+                       CALL flioinqa( knumfl, cdvar, 'offset', ll_fnd )
+                       IF( ll_fnd ) THEN
+                          CALL fliogeta( knumfl, cdvar, 'offset', iom_file(knumfl)%ofs(idvar) )
+                       ELSE
+                          iom_file(knumfl)%ofs(idvar) = 0.
+                       END IF
+                       !
+                       iom_varid = idvar
+                       IF ( PRESENT(kdimsz) ) THEN 
+                          IF ( i_nvd == size(kdimsz) ) THEN
+                             kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar)
+                          ELSE
+                             WRITE(ctmp1,*) i_nvd, size(kdimsz)
+                             CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) )
+                          ENDIF
+                       ENDIF
+                    ELSE
+                       CALL ctl_stop( trim(clinfo), 'Too many dimensions in the file '//iom_file(knumfl)%name,   &
+                            &           'increase the parameter jpmax_vars')
+                    ENDIF
+                 ELSE  
+                    CALL ctl_warn( trim(clinfo), 'Variable '//trim(cdvar)// &
+                         & ' is not found in the file '//trim(iom_file(knumfl)%name) )
+                 ENDIF
+              ELSE
+                 CALL ctl_stop( trim(clinfo), 'Too many variables in the file '//iom_file(knumfl)%name,   &
+                      &           'increase the parameter jpmax_vars')
+              ENDIF
+           ELSE
+              iom_varid = idvar
+              IF ( PRESENT(kdimsz) ) THEN 
+                 i_nvd = iom_file(knumfl)%ndims(idvar)
+                 IF ( i_nvd == size(kdimsz) ) THEN
+                    kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar)
+                 ELSE
+                    WRITE(ctmp1,*) i_nvd, size(kdimsz)
+                    CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) )
+                 ENDIF
+              ENDIF
+           ENDIF
+        ENDIF
+     ENDIF
+
+   END FUNCTION iom_varid
+
+  
+   FUNCTION iom_get_gblatt( knumfl, kinfonum )
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  iom_get_gblatt  ***
+      !!
+      !! ** Purpose : ???
+      !!
+      !! ** Method : ???
+      !!
+      !!-----------------------------------------------------------------------
+      INTEGER,INTENT(in) :: knumfl
+      INTEGER,          intent(in) :: kinfonum
+      !
+      CHARACTER(LEN=10) :: clinfo
+      REAL(wp) :: iom_get_gblatt
+      !!-----------------------------------------------------------------------
+
+      WRITE(clinfo,*) kinfonum
+      clinfo = 'info'//trim(adjustl(clinfo))
+      CALL fliogeta (knumfl, "?", clinfo, iom_get_gblatt)
+
+   END FUNCTION iom_get_gblatt
+
+   !!======================================================================
+END MODULE iom
