Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 5602)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 6772)
@@ -44,9 +44,10 @@
    USE sbc_ice         ! Surface boundary condition: ice fields
    USE lib_fortran     ! to use key_nosignedzero
+   USE sbcapr
 #if defined key_lim3
    USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b
    USE limthd_dh       ! for CALL lim_thd_snwblow
 #elif defined key_lim2
-   USE ice_2, ONLY     : u_ice, v_ice
+   USE ice_2, ONLY     : u_ice, v_ice, pfrld
    USE par_ice_2
 #endif
@@ -83,4 +84,6 @@
    REAL(wp), PARAMETER ::   Cice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice
    REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant
+   REAL(wp), PARAMETER ::   rgas =  287.1         ! gas const. dry air (J/kg/K)
+   REAL(wp), PARAMETER ::   rvap =  461.51        ! gas const. vapour  (J/kg/K)
 
    !                                  !!* Namelist namsbc_core : CORE bulk parameters
@@ -91,4 +94,11 @@
    REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
+   !
+   LOGICAL  ::   ln_tair_celsius  !: logical flag for Read Tair: Tair in NEMO is Kelvin
+   LOGICAL  ::   ln_humi_rel      !: logical flag for Read relative humidity (T) or specific humidity (F)
+   LOGICAL  ::   ln_cohum_arc     !: logical flag for Correction of Humidity in the Arctic Ocean
+   LOGICAL  ::   ln_cotair_arc    !: logical flag for Correction of Air Temperature in the Arctic Ocean
+   LOGICAL  ::   ln_corad_antar   !: logical flag for Correction of radiatives fluxes in the Southern Ocean
+
 
    !! * Substitutions
@@ -143,4 +153,10 @@
       INTEGER  ::   ios      ! Local integer output status for namelist read
       !
+      INTEGER  ::   ji,jj
+      REAL(wp) ::   zzlat, zzlat1, zzlat2, zfm, zfrld
+      REAL(wp) ::   zmin,zmax
+      REAL(wp), DIMENSION(:,:), POINTER :: xyt,z_qsr,z_qlw,z_qsr1,z_qlw1, z_hum, z_tair
+      REAL(wp), DIMENSION(:,:), POINTER :: zqsr_lr, zqsr_hr, zqlw_lr, zqlw_hr
+     
       CHARACTER(len=100) ::  cn_dir   !   Root directory for location of core files
       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
@@ -151,7 +167,12 @@
          &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
          &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
-         &                  sn_tdif, rn_zqt,  rn_zu
+         &                  sn_tdif, rn_zqt,  rn_zu , ln_tair_celsius,   &
+         &                  ln_humi_rel  , ln_cohum_arc,      &
+         &                  ln_cotair_arc, ln_corad_antar
+
       !!---------------------------------------------------------------------
       !
+      CALL wrk_alloc( jpi,jpj, xyt,z_qsr,z_qlw,z_qsr1,z_qlw1, z_hum, z_tair )
+      CALL wrk_alloc( jpi,jpj, zqsr_lr, zqsr_hr, zqlw_lr, zqlw_hr )
       !                                         ! ====================== !
       IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
@@ -194,4 +215,13 @@
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
          !
+         !
+         IF(lwp) WRITE(numout,*) 'sbc_blk_core: jfld = ',jfld
+         IF( ln_cohum_arc   ) CALL ctl_warn( 'sbc_blk_core: correction of humidity in arctic' )
+         IF( ln_cotair_arc  ) CALL ctl_warn( 'sbc_blk_core: correction of air temperature in arctic' )
+         IF( ln_corad_antar ) CALL ctl_warn( 'sbc_blk_core: correction of short radiation in antartic' )
+         IF( ln_humi_rel    ) CALL ctl_warn( 'sbc_blk_core: use relative humidity instead of specific humidity')
+         IF( ln_tair_celsius) CALL ctl_warn( 'sbc_blk_core: Tair is read in Celsius')
+         IF(lwp) WRITE(numout,*) 'sbc_blk_core: rn_pfac = ',rn_pfac
+         !
          sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
          !
@@ -199,4 +229,102 @@
 
       CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
+
+      !=========================================
+      !  ONLINE CORRECTIONS
+      !=========================================
+      !
+      ! Correction of Tair
+      !
+      IF( ln_tair_celsius .AND. MOD( kt-1, nn_fsbc ) == 0 ) THEN
+         sf(jp_tair)%fnow = sf(jp_tair)%fnow + 273.15_wp  ! Conversion of the Temperature °C --> Kelvin
+      ENDIF
+      !
+      ! Correction of SW and LW in the Southern Ocean
+      !
+      IF( ln_corad_antar .AND. .NOT. sf(jp_qsr)%ln_tint .AND. MOD( kt-1, 86400/INT(rdt) ) == 0 ) THEN
+         z_qsr(:,:) = 0.8 * sf(jp_qsr)%fnow(:,:,1)
+         xyt(:,:) = 0.e0 ; zzlat1 = -65. ; zzlat2 = -60.
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zzlat = gphit(ji,jj)
+               IF( zzlat >= zzlat1 .AND. zzlat <= zzlat2 ) THEN
+                  xyt(ji,jj) = (zzlat2-zzlat)/(zzlat2-zzlat1)
+               ELSE IF ( zzlat < zzlat1 ) THEN
+                  xyt(ji,jj) = 1
+               ENDIF
+            END DO
+         END DO
+         IF(lwp) WRITE(numout,*) 'Correc ln_corad_antar'
+         z_qsr1(:,:) = z_qsr(:,:) * xyt(:,:) + ( 1.0 - xyt(:,:) ) * sf(jp_qsr)%fnow(:,:,1)
+         sf(jp_qsr)%fnow(:,:,1) = z_qsr1(:,:)
+      ENDIF
+
+      IF( MOD( kt-1, nn_fsbc ) == 0 )THEN
+         !
+         IF ( nmonth >= 5 .AND. nmonth <= 9 ) THEN
+            !
+            ! Correction of Humidity in the Arctic Ocean
+            !
+            IF( ln_cohum_arc ) THEN
+               z_hum(:,:) = 0.85 * sf(jp_humi)%fnow(:,:,1)
+               xyt(:,:) = 0.e0 ; zzlat1 = 78. ; zzlat2 = 82.
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zzlat = gphit(ji,jj)
+#if defined key_lim2 ||  defined key_lim3 
+                     IF ( ALLOCATED(pfrld) ) THEN ; zfrld = pfrld(ji,jj) ; ELSE ; zfrld = 0 ; ENDIF
+#endif
+                     IF( zzlat >= zzlat1 .AND. zzlat <= zzlat2 .AND. zfrld < 0.85 ) THEN
+                        xyt(ji,jj) = ( zzlat - zzlat1 ) / ( zzlat2 - zzlat1 )
+                     ELSE IF ( zzlat > zzlat2 .AND. zfrld < 0.85 ) THEN
+                        xyt(ji,jj) = 1._wp
+                     ENDIF
+                  ENDDO
+               ENDDO
+               IF(lwp) WRITE(numout,*) 'Correc ln_cohum_arc'
+               sf(jp_humi)%fnow(:,:,1) = z_hum(:,:) * xyt(:,:) + ( 1.0 - xyt(:,:) ) * sf(jp_humi)%fnow(:,:,1)
+            ENDIF
+            !
+            ! Correction of Air Temperature in the Arctic Ocean
+            !
+            IF( ln_cotair_arc ) THEN
+               z_tair(:,:) = sf(jp_tair)%fnow(:,:,1) - 2.0
+               xyt(:,:) = 0.e0 ; zzlat1 = 78. ; zzlat2 = 82.
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zzlat = gphit(ji,jj)
+#if defined key_lim2 ||  defined key_lim3 
+                     IF( ALLOCATED(pfrld) ) THEN ; zfrld = pfrld(ji,jj) ; ELSE ; zfrld=0 ; ENDIF
+#endif
+                     IF( zzlat >= zzlat1 .AND. zzlat <= zzlat2 .AND. zfrld < 0.85 ) THEN
+                        xyt(ji,jj) = ( zzlat - zzlat1 ) / ( zzlat2 - zzlat1 )
+                     ELSE IF( zzlat > zzlat2 .AND. zfrld < 0.85 ) THEN
+                        xyt(ji,jj) = 1._wp
+                     ENDIF
+                  END DO
+               ENDDO
+               IF(lwp) WRITE(numout,*) 'Correc ln_cotair_arc'
+               sf(jp_tair)%fnow(:,:,1) = z_tair(:,:) * xyt(:,:) + ( 1.0 - xyt(:,:) ) * sf(jp_tair)%fnow(:,:,1)
+            ENDIF
+            !
+         ENDIF ! 5 <= nmonth <= 9
+
+         !
+      ENDIF ! IF MOD( kt-1, nn_fsbc )
+
+      DO jj=1,jpj
+         DO ji=1,jpi
+            sf(jp_humi)%fnow(ji,jj,1) = MAX( MIN( sf(jp_humi)%fnow(ji,jj,1) ,1.0 ) , 0.0 )
+            sf(jp_prec)%fnow(ji,jj,1) = MAX(      sf(jp_prec)%fnow(ji,jj,1) ,0.0 )
+            sf(jp_qsr )%fnow(ji,jj,1) = MAX(      sf(jp_qsr )%fnow(ji,jj,1) ,0.0 )
+            sf(jp_qlw )%fnow(ji,jj,1) = MAX(      sf(jp_qlw )%fnow(ji,jj,1) ,0.0 )
+         ENDDO
+      END DO
+
+      !
+      !=========================================
+      ! END OF ONLINE CORRECTIONS
+      !=========================================
+      !
 
       !                                            ! compute the surface ocean fluxes using CORE bulk formulea
@@ -215,4 +343,7 @@
       ENDIF
 #endif
+      !
+      CALL wrk_dealloc( jpi,jpj, xyt,z_qsr,z_qlw,z_qsr1,z_qlw1, z_hum, z_tair )
+      CALL wrk_dealloc( jpi,jpj, zqsr_lr, zqsr_hr, zqlw_lr, zqlw_hr )
       !
    END SUBROUTINE sbc_blk_core
@@ -257,4 +388,6 @@
       REAL(wp), DIMENSION(:,:), POINTER ::   zt_zu             ! air temperature at wind speed height
       REAL(wp), DIMENSION(:,:), POINTER ::   zq_zu             ! air spec. hum.  at wind speed height
+      REAL(wp), DIMENSION(:,:), POINTER ::   zqatm , zpatm     ! specific humidity and mean sea level pressure (Pa)
+      REAL(wp) :: vt, vp, vq, zqa, zq0, zq1, zq2, zee
       !!---------------------------------------------------------------------
       !
@@ -262,5 +395,5 @@
       !
       CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap )
-      CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
+      CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu ,zqatm, zpatm )
       !
       ! local scalars ( place there for vector optimisation purposes)
@@ -314,7 +447,30 @@
       ! ... specific humidity at SST and IST
       zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) )
-
+      !
+      IF ( ln_humi_rel ) THEN
+         zq0    = rvap / rgas - 1.0
+         zq1    = rgas / rvap
+         zq2    = 1.0 - zq1
+         zpatm(:,:) = 100800.                   ! atmospheric pressure (assumed constant  here)
+         IF ( ln_apr_dyn ) zpatm(:,:) = apr(:,:)
+         DO jj = 1 , jpj
+            DO ji = 1 , jpi
+               vt  = sf(jp_tair)%fnow(ji,jj,1) - rt0  ! air temperature (Celsius)
+               vp  = zpatm(ji,jj) / 100.              ! mean sea level pressure (mb or hPa)
+               vq  = sf(jp_humi)%fnow(ji,jj,1)        ! relative humidity (fraction of 1)
+               ! Convert RH at the air/sea interface in specific humidity (kg/kg)
+               ! Teten's formula for qsat (mb)
+               zqa = ( 1.0007 + 3.46e-6 * vp) * 6.1121 * EXP( 17.502 * vt / ( 240.97+vt ) )
+               zee = zqa * vq                         ! vapour partial pressure (mb)
+               vq  = zq1 * zee / ( vp - zq2 * zee )   ! specific humidity (kg/kg)
+               zqatm(ji,jj) = vq
+            ENDDO
+         ENDDO
+      ELSE
+         zqatm(:,:)=sf(jp_humi)%fnow(:,:,1)
+      ENDIF
+      !
       ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point :
-      CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   &
+      CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, zqatm, wndm,   &
          &               Cd, Ch, Ce, zt_zu, zq_zu )
     
@@ -354,5 +510,6 @@
       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
          !! q_air and t_air are (or "are almost") given at 10m (wind reference height)
-         zevap(:,:) = rn_efac*MAX( 0._wp,     rhoa*Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) )*wndm(:,:) ) ! Evaporation
+         !zevap(:,:) = rn_efac*MAX( 0._wp,     rhoa*Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) )*wndm(:,:) ) ! Evaporation
+          zevap(:,:) = rn_efac*MAX( 0._wp,     rhoa*Ce(:,:)*( zqsatw(:,:) - zqatm(:,:)              )*wndm(:,:) ) ! Evaporation
          zqsb (:,:) =                     cpa*rhoa*Ch(:,:)*( zst   (:,:) - sf(jp_tair)%fnow(:,:,1) )*wndm(:,:)   ! Sensible Heat
       ELSE
@@ -414,5 +571,5 @@
       !
       CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap )
-      CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
+      CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu, zqatm, zpatm )
       !
       IF( nn_timing == 1 )  CALL timing_stop('blk_oce_core')
@@ -437,5 +594,5 @@
       REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point
       !!---------------------------------------------------------------------
-      !
+
       IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau')
       !
@@ -530,4 +687,5 @@
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla, zcoef_dqsb
       REAL(wp) ::   zztmp, z1_lsub                               ! temporary variable
+      REAL(wp) ::   ztamr,zmt1,zmt2,zmt3,zev,zes
       !!
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice
@@ -536,4 +694,5 @@
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw       ! evaporation and snw distribution after wind blowing (LIM3)
+      REAL(wp), DIMENSION(:,:)  , POINTER ::   zqatm, zpatm , ztatm            ! specific humidity
       !!---------------------------------------------------------------------
       !
@@ -541,4 +700,25 @@
       !
       CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
+      CALL wrk_alloc( jpi,jpj, zqatm, zpatm, ztatm )
+ 
+     IF ( ln_humi_rel ) THEN
+         zpatm(:,:) = 100800.                   ! atmospheric pressure (assumed constant here)
+         IF (ln_apr_dyn) zpatm(:,:) = apr(:,:)
+         DO jj=1,jpj
+            DO ji=1,jpi
+               ztatm (ji,jj) = sf(jp_tair)%fnow(ji,jj,1)                   ! air temperature in Kelvins
+               ztamr = ztatm(ji,jj) - rtt                                  ! Saturation water vapour
+               zmt1  = SIGN( 17.269,  ztamr )
+               zmt2  = SIGN( 21.875,  ztamr )
+               zmt3  = SIGN( 28.200, -ztamr )
+               zes   = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &
+                  &                / ( ztatm(ji,jj) - 35.86  + MAX( 0.e0, zmt3 ) )  )
+               zev = sf(jp_humi)%fnow(ji,jj,1) * zes                       ! vapour pressure
+               zqatm(ji,jj) = 0.622 * zev / ( zpatm(ji,jj) - 0.378 * zev ) ! specific humidity
+            ENDDO
+         ENDDO
+      ELSE
+         zqatm(:,:) = sf(jp_humi)%fnow(:,:,1)
+      ENDIF
 
       ! local scalars ( place there for vector optimisation purposes)
@@ -574,5 +754,5 @@
                ! Latent Heat
                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                           
-                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
+                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - zqatm(ji,jj)  ) )
               ! Latent heat sensitivity for ice (Dqla/Dt)
                IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
@@ -659,4 +839,5 @@
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx')
       
+      CALL wrk_dealloc( jpi,jpj, zqatm, zpatm, ztatm )
    END SUBROUTINE blk_ice_core_flx
 #endif
