Index: branches/2011/dev_r2784_CMCC1_topbfm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/2011/dev_r2784_CMCC1_topbfm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 2798)
+++ branches/2011/dev_r2784_CMCC1_topbfm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 2799)
@@ -44,5 +44,4 @@
    PUBLIC   blk_ice_core         ! routine called in sbc_ice_lim module
 
-   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read 
    INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
    INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
@@ -54,4 +53,12 @@
    INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
    INTEGER , PARAMETER ::   jp_tdif = 9           ! index of tau diff associated to HF tau   (N/m2)   at T-point
+#if defined key_orca_r025
+   INTEGER , PARAMETER ::   jp_swc  = 10          ! index of GEWEX correction for SW radiation  at T-point
+   INTEGER , PARAMETER ::   jp_lwc  = 11          ! index of GEWEX correction for LW radiation  at T-point
+   INTEGER , PARAMETER ::   jp_prc  = 12          ! index of PMWC correction forat T-point
+   INTEGER , PARAMETER ::   jpfld   = 12          ! maximum number of files to read
+#else
+   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read
+#endif
    
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
@@ -71,4 +78,23 @@
    REAL(wp) ::   rn_pfac   = 1.        ! multiplication factor for precipitation
 
+#if defined key_orca_r025
+   LOGICAL  ::   ln_printdia= .TRUE.     ! logical flag for height of air temp. and hum
+   LOGICAL  ::   ln_netsw   = .TRUE.     ! logical flag for height of air temp. and hum
+   LOGICAL  ::   ln_core_graceopt=.FALSE., ln_core_spinup=.FALSE.
+   LOGICAL  ::   ln_gwxc = .TRUE.
+   LOGICAL  ::   ln_corad_antar =.FALSE., ln_corad_arc =.FALSE. , ln_cotair_arc = .FALSE.
+   LOGICAL  ::   ln_coprecip =.FALSE.
+   REAL(wp) ::   rn_qns_bias = 0._wp     ! heat flux bias
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  area 
+   REAL(wp) :: araux
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  zqlw, zqsb  ! long wave and sensible heat fluxes
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  zqla, zevap ! latent heat fluxes and evaporation
+
+   REAL(wp), PARAMETER :: zalph = 2.408724e-06_wp, &
+   & zbet = -0.006936579_wp, zgam = 449.9094_wp ! GRACE regression coefficients
+#endif
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -80,5 +106,5 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
+   
    SUBROUTINE sbc_blk_core( kt )
       !!---------------------------------------------------------------------
@@ -112,4 +138,7 @@
       !!              - emp, emps   evaporation minus precipitation
       !!----------------------------------------------------------------------
+#if defined key_orca_r025
+      USE ice_2
+#endif
       INTEGER, INTENT(in) ::   kt   ! ocean time step
       !!
@@ -117,4 +146,5 @@
       INTEGER  ::   ifpr     ! dummy loop indice
       INTEGER  ::   jfld     ! dummy loop arguments
+      INTEGER  ::   ji, jj
       !!
       CHARACTER(len=100) ::  cn_dir   !   Root directory for location of core files
@@ -123,12 +153,34 @@
       TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !   "                                 "
       TYPE(FLD_N) ::   sn_tdif                                 !   "                                 "
+      TYPE(FLD_N) ::   sn_swc, sn_lwc                          !   "                                 "
+      TYPE(FLD_N) ::   sn_prc
+      INTEGER  ::   iter_shapiro = 250
+      REAL :: zzlat, zzlat1, zzlat2, zfm, zfrld, ztmp
+      REAL(wp), DIMENSION(jpi,jpj):: xyt,z_qsr,z_qlw,z_qsr1,z_qlw1, z_hum, z_tair
+      REAL(wp), DIMENSION(jpi,jpj):: zqsr_lr, zqsr_hr, zqlw_lr, zqlw_hr, zprec_hr, zprec_lr
+      CHARACTER(len=20)  ::  c_kind='ORCA_GLOB'
+#if defined key_orca_r025
+      NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac,           &
+         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
+         &                  sn_qlw , sn_tair, sn_prec  , sn_snow, sn_tdif   &
+         &                  sn_swc , sn_lwc , sn_prc   , ln_gwxc  ,         &
+         &                  ln_corad_antar, ln_corad_arc, ln_cotair_arc, ln_coprecip ,  &
+         &                  rn_qns_bias, ln_printdia, ln_netsw, ln_core_graceopt,ln_core_spinup
+      !!---------------------------------------------------------------------
+#else
       NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac,           &
          &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
          &                  sn_qlw , sn_tair, sn_prec  , sn_snow, sn_tdif
       !!---------------------------------------------------------------------
+#endif
 
       !                                         ! ====================== !
       IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
          !                                      ! ====================== !
+#if defined key_orca_r025
+         !                                      !==  allocate sbc arrays
+            IF( sbc_blk_core_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_core_alloc : unable to allocate arrays' )
+#endif
+
          ! set file information (default values)
          cn_dir = './'       ! directory in which the model is executed
@@ -146,4 +198,9 @@
          sn_snow = FLD_N( 'snow'   ,    -1     , 'snow'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       )
          sn_tdif = FLD_N( 'taudif' ,    24     , 'taudif' ,  .true.    , .false. ,   'yearly'  , ''       , ''       )
+#if defined key_orca_r025
+         sn_swc  = FLD_N( 'swc'    ,    24     ,  'swc'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       )
+         sn_lwc  = FLD_N( 'lwc'    ,    24     ,  'lwc'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       )
+         sn_prc  = FLD_N( 'prc'    ,    24     ,  'prc'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       )
+#endif
          !
          REWIND( numnam )                          ! read in namlist namsbc_core
@@ -163,7 +220,16 @@
          slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
          slf_i(jp_tdif) = sn_tdif
+#if defined key_orca_r025
+            slf_i(jp_swc ) = sn_swc
+            slf_i(jp_lwc ) = sn_lwc
+            slf_i(jp_prc ) = sn_prc
+#endif
          !                 
          lhftau = ln_taudif                        ! do we use HF tau information?
          jfld = jpfld - COUNT( (/.NOT. lhftau/) )
+#if defined key_orca_r025
+         IF( .NOT. ln_gwxc )     jfld = jfld - 2
+         IF( .NOT. ln_coprecip ) jfld = jfld - 1
+#endif
          !
          ALLOCATE( sf(jfld), STAT=ierror )         ! set sf structure
@@ -176,11 +242,125 @@
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
          !
+#if defined key_orca_r025
+         IF( ln_printdia .OR. ln_core_graceopt ) THEN
+           area = (e1t * e2t)
+           araux = sum ( area * tmask(:,:,1) )
+           IF(lk_mpp) CALL mpp_sum ( araux )
+         ENDIF
+#endif
       ENDIF
 
       CALL fld_read( kt, nn_fsbc, sf )        ! input fields provided at the current time-step
 
-      !                                                        ! surface ocean fluxes computed with CLIO bulk formulea
-      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m )
-      !
+      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
+
+#if defined key_orca_r025
+      ! Introduce ERA-Interim filtering and correction
+
+         IF( ln_gwxc ) THEN
+
+           call Shapiro_1D(sf(jp_qsr)%fnow(:,:,1),iter_shapiro, c_kind, zqsr_lr)
+           zqsr_hr(:,:)=sf(jp_qsr)%fnow(:,:,1)-zqsr_lr(:,:)          ! We get large scale and small scale
+
+           call Shapiro_1D(sf(jp_qlw)%fnow(:,:,1),iter_shapiro, c_kind, zqlw_lr)
+           zqlw_hr(:,:)=sf(jp_qlw)%fnow(:,:,1)-zqlw_lr(:,:)          ! We get large scale and small scale
+
+           z_qsr1(:,:)=zqsr_lr(:,:)*sf(jp_swc)%fnow(:,:,1) + zqsr_hr(:,:)
+           z_qlw1(:,:)=zqlw_lr(:,:)*sf(jp_lwc)%fnow(:,:,1) + zqlw_hr(:,:)
+
+           DO jj=1,jpj
+             DO ji=1,jpi
+               z_qsr1(ji,jj)=max(z_qsr1(ji,jj),0.0)
+               z_qlw1(ji,jj)=max(z_qlw1(ji,jj),0.0)
+             END DO
+           END DO
+
+         ENDIF
+
+         IF( ln_coprecip ) THEN
+
+           call Shapiro_1D(sf(jp_prec)%fnow(:,:,1),iter_shapiro,c_kind,zprec_lr)
+           zprec_hr(:,:)=sf(jp_prec)%fnow(:,:,1)-zprec_lr(:,:)       ! We get large scale and small scale
+
+           DO jj=1,jpj
+             DO ji=1,jpi
+               IF( zprec_lr(ji,jj) .GT. 0._wp ) THEN
+                  ztmp = LOG( ( 1000._wp + sf(jp_prc)%fnow(ji,jj,1) ) * EXP( zprec_lr(ji,jj) ) / 1000._wp )
+                  sf(jp_prec)%fnow(ji,jj,1) = max(ztmp+zprec_hr(ji,jj),0.0)
+               ENDIF
+             END DO
+           END DO
+
+         ENDIF
+
+         IF ( ln_corad_antar ) THEN           ! correction of SW and LW in the Southern Ocean
+ 
+           z_qsr(:,:)=0.8*z_qsr1(:,:)
+           z_qlw(:,:)=1.1*z_qlw1(:,:)
+           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
+           z_qsr1(:,:)=z_qsr(:,:)*xyt(:,:)+(1.0-xyt(:,:))*z_qsr1(:,:)
+           z_qlw1(:,:)=z_qlw(:,:)*xyt(:,:)+(1.0-xyt(:,:))*z_qlw1(:,:)
+
+         ENDIF
+
+         IF ( ln_corad_arc ) THEN         ! correction of SW in the Arctic Ocean
+
+           z_qsr(:,:)=0.7*z_qsr1(:,:)
+           xyt(:,:) = 0.e0
+           zzlat1 = 78.
+           zzlat2 = 82.
+           DO jj = 1, jpj
+             DO ji = 1, jpi
+               zzlat = gphit(ji,jj)
+               IF ( zzlat >= zzlat1 .AND. zzlat <= zzlat2 ) THEN
+                  xyt(ji,jj) = (zzlat-zzlat1)/(zzlat2-zzlat1)
+               ELSE IF ( zzlat > zzlat2 ) THEN
+                  xyt(ji,jj) = 1
+               ENDIF
+             END DO
+           END DO
+           z_qsr1(:,:)=z_qsr(:,:)*xyt(:,:)+(1.0-xyt(:,:))*z_qsr1(:,:)
+
+         ENDIF
+
+         sf(jp_qsr)%fnow(:,:,1)=z_qsr1(:,:)
+         sf(jp_qlw)%fnow(:,:,1)=z_qlw1(:,:)
+
+         IF ( ln_cotair_arc ) THEN           ! correction of Air Temperature in the Arctic Ocean
+
+           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) ; zfrld=frld(ji,jj)
+               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
+               ENDIF
+             END DO
+           END DO
+           sf(jp_tair)%fnow(:,:,1)=z_tair(:,:)*xyt(:,:)+(1.0-xyt(:,:))*sf(jp_tair)%fnow(:,:,1)
+
+         ENDIF
+
+#endif
+
+         CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m )   ! compute the surface ocean fluxes using CLIO bulk formulea
+
+      ENDIF
+      !                                                  ! using CORE bulk formulea
    END SUBROUTINE sbc_blk_core
    
@@ -315,4 +495,9 @@
       IF( lhftau ) THEN 
 !CDIR COLLAPSE
+#if defined key_orca_r025
+         ! Changed!!! Multiply by QSCAT correction
+           zwnd_i(:,:) = zwnd_i(:,:) * sf(jp_tdif)%fnow(:,:,1)
+           zwnd_j(:,:) = zwnd_j(:,:) * sf(jp_tdif)%fnow(:,:,1)
+#endif
          taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
       ENDIF
@@ -938,4 +1123,104 @@
       !
     END FUNCTION psi_h
+
+#if defined key_orca_r025
+    INTEGER FUNCTION sbc_blk_core_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE sbc_blk_core_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( area(jpi,jpj)         , zqlw(jpi,jpj)      ,     &
+         &      zqsb(jpi,jpj)         , zqla(jpi,jpj)      ,     &
+         &      zevap(jpi,jpj)        , STAT=sbc_blk_core_alloc )
+      zqlw=0._wp
+      zqsb=0._wp
+	!
+      IF( lk_mpp            )   CALL mpp_sum ( sbc_blk_core_alloc )
+      IF( sbc_blk_core_alloc > 0 )   CALL ctl_warn('sbc_blk_core_alloc: allocation of arrays failed')
+    END FUNCTION sbc_blk_core_alloc
+#endif
+
+    SUBROUTINE Shapiro_1D(rla_varin,id_np, cd_overlap, rlpa_varout) !GIG
+      !!=====================================================================
+      !!
+      !! Description: This function applies a 1D Shapiro filter
+      !!              (3 points filter) horizontally to a 2D field
+      !!              in regular grid
+      !! Arguments :
+      !!            rla_varin   : Input variable to filter
+      !!            zla_mask    : Input mask variable
+      !!            id_np       : Number of Shapiro filter iterations
+      !!            cd_overlap  : Logical argument for periodical condition
+      !!                          (global ocean case)
+      !!            rlpa_varout : Output filtered variable
+      !!
+      !! History : 08/2009  S. CAILLEAU : from 1st version of N. FERRY
+      !!           09/2009  C. REGNIER  : Corrections
+      !!
+      !!=====================================================================
+      IMPLICIT NONE
+      INTEGER, INTENT(IN)                       :: id_np
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)  :: rla_varin !GIG
+      CHARACTER(len=20), INTENT(IN)             :: cd_overlap !GIG
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT) :: rlpa_varout !GIG
+
+      REAL(wp), DIMENSION(jpi,jpj)              :: rlpa_varout_tmp
+      REAL, PARAMETER                           :: rl_alpha = 1./2.    ! fixed stability coefficient (isotrope case)
+      REAL, parameter                           :: rap_aniso_diff_XY=2.25 ! anisotrope case
+      REAL                                      :: alphax,alphay, znum, zden,test
+      INTEGER                                   :: ji, jj, jn, nn
+!
+!! rap_aniso_diff_XY=2.25 : valeur trouvée empiriquement pour 140 itération pour le filtre de Shapiro et
+!! pour un rapport d'anisotopie de 1.5 : on filtre de plus rapidement en x qu'eny.
+!------------------------------------------------------------------------------
+!
+! Loop on several filter iterations
+
+!     Global ocean case
+      IF (( cd_overlap == 'MERCA_GLOB' )   .OR.   &
+          ( cd_overlap == 'REGULAR_GLOB' ) .OR.   &
+          ( cd_overlap == 'ORCA_GLOB' )) THEN
+             rlpa_varout(:,:) = rla_varin(:,:)
+             rlpa_varout_tmp(:,:) = rlpa_varout(:,:)
+!
+
+       alphax=1./2.
+       alphay=1./2.
+!  Dx/Dy=rap_aniso_diff_XY  , D_ = vitesse de diffusion
+!  140 passes du fitre, Lx/Ly=1.5, le rap_aniso_diff_XY correspondant est:
+       IF ( rap_aniso_diff_XY .GE. 1. ) alphay=alphay/rap_aniso_diff_XY
+       IF ( rap_aniso_diff_XY .LT. 1. ) alphax=alphax*rap_aniso_diff_XY
+
+        DO jn = 1,id_np   ! number of passes of the filter
+            DO ji = 2,jpim1
+               DO jj = 2,jpjm1
+                  ! We crop on the coast
+                   znum = rlpa_varout_tmp(ji,jj)   &
+                          + 0.25*alphax*(rlpa_varout_tmp(ji-1,jj  )-rlpa_varout_tmp(ji,jj))*tmask(ji-1,jj  ,1)  &
+                          + 0.25*alphax*(rlpa_varout_tmp(ji+1,jj  )-rlpa_varout_tmp(ji,jj))*tmask(ji+1,jj  ,1)  &
+                          + 0.25*alphay*(rlpa_varout_tmp(ji  ,jj-1)-rlpa_varout_tmp(ji,jj))*tmask(ji  ,jj-1,1)  &
+                          + 0.25*alphay*(rlpa_varout_tmp(ji  ,jj+1)-rlpa_varout_tmp(ji,jj))*tmask(ji  ,jj+1,1)
+                   rlpa_varout(ji,jj)=znum*tmask(ji,jj,1)+rla_varin(ji,jj)*(1.-tmask(ji,jj,1))
+                ENDDO  ! end loop ji
+            ENDDO  ! end loop jj
+!
+!
+!           Periodical condition in case of cd_overlap (global ocean)
+!           - on a mercator projection grid we consider that singular point at poles
+!             are a mean of the values at points of the previous latitude
+!           - on ORCA and regular grid we copy the values at points of the previous latitude
+            IF ( cd_overlap == 'MERCAT_GLOB' ) THEN
+!GIG case unchecked
+               rlpa_varout(1,1) = SUM(rlpa_varout(:,2)) / jpi
+               rlpa_varout(jpi,jpj) = SUM(rlpa_varout(:,jpj-1)) / jpi
+            ELSE
+               call lbc_lnk(rlpa_varout, 'T', 1.) ! Boundary condition
+            ENDIF
+            rlpa_varout_tmp(:,:) = rlpa_varout(:,:)
+         ENDDO  ! end loop jn
+      ENDIF
+
+!
+    END SUBROUTINE Shapiro_1D
+
   
    !!======================================================================
