Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 4627)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 4990)
@@ -17,4 +17,5 @@
    USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic
    USE in_out_manager  ! I/O manager
+   USE iom, ONLY : iom_put,iom_use              ! I/O manager library !!Joakim edit
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -23,5 +24,4 @@
    USE daymod          ! calendar
    USE fldread         ! read input fields
-
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice   fields
@@ -38,4 +38,5 @@
    USE ice_calendar, only: dt
    USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
+# if defined key_cice4
    USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
                 sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm,     &
@@ -44,7 +45,16 @@
                 uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
                 swvdr,swvdf,swidr,swidf
+   USE ice_therm_vertical, only: calc_Tsfc
+#else
+   USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
+                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,     &
+                fresh_ai,fhocn_ai,fswthru_ai,frzmlt,          &
+                flatn_f,fsurfn_f,fcondtopn_f,                    &
+                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
+                swvdr,swvdf,swidr,swidf
+   USE ice_therm_shared, only: calc_Tsfc
+#endif
    USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
    USE ice_atmo, only: calc_strair
-   USE ice_therm_vertical, only: calc_Tsfc
 
    USE CICE_InitMod
@@ -95,5 +105,5 @@
    END FUNCTION sbc_ice_cice_alloc
 
-   SUBROUTINE sbc_ice_cice( kt, nsbc )
+   SUBROUTINE sbc_ice_cice( kt, ksbc )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_ice_cice  ***
@@ -113,5 +123,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt      ! ocean time step
-      INTEGER, INTENT(in) ::   nsbc    ! surface forcing type
+      INTEGER, INTENT(in) ::   ksbc    ! surface forcing type
       !!----------------------------------------------------------------------
       !
@@ -123,15 +133,15 @@
 
          ! Make sure any fluxes required for CICE are set
-         IF ( nsbc == 2 )  THEN
+         IF      ( ksbc == jp_flx ) THEN
             CALL cice_sbc_force(kt)
-         ELSE IF ( nsbc == 5 ) THEN
+         ELSE IF ( ksbc == jp_cpl ) THEN
             CALL sbc_cpl_ice_flx( 1.0-fr_i  )
          ENDIF
 
-         CALL cice_sbc_in ( kt, nsbc )
+         CALL cice_sbc_in  ( kt, ksbc )
          CALL CICE_Run
-         CALL cice_sbc_out ( kt, nsbc )
-
-         IF ( nsbc == 5 )  CALL cice_sbc_hadgam(kt+1)
+         CALL cice_sbc_out ( kt, ksbc )
+
+         IF ( ksbc == jp_cpl )  CALL cice_sbc_hadgam(kt+1)
 
       ENDIF                                          ! End sea-ice time step only
@@ -141,13 +151,13 @@
    END SUBROUTINE sbc_ice_cice
 
-   SUBROUTINE cice_sbc_init (nsbc)
+   SUBROUTINE cice_sbc_init (ksbc)
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE cice_sbc_init  ***
       !! ** Purpose: Initialise ice related fields for NEMO and coupling
       !!
-      INTEGER, INTENT( in  ) ::   nsbc                ! surface forcing type
+      INTEGER, INTENT( in  ) ::   ksbc                ! surface forcing type
       REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
       REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar
-      INTEGER  ::   ji, jj, jl                        ! dummy loop indices
+      INTEGER  ::   ji, jj, jl, jk                    ! dummy loop indices
       !!---------------------------------------------------------------------
 
@@ -161,13 +171,22 @@
       jj_off = INT ( (jpjglo - ny_global) / 2 )
 
+#if defined key_nemocice_decomp
+      ! Pass initial SST from NEMO to CICE so ice is initialised correctly if
+      ! there is no restart file.
+      ! Values from a CICE restart file would overwrite this
+      IF ( .NOT. ln_rstart ) THEN    
+         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
+      ENDIF  
+#endif
+
 ! Initialize CICE
       CALL CICE_Initialize
 
 ! Do some CICE consistency checks
-      IF ( (nsbc == 2) .OR. (nsbc == 5) ) THEN
+      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_cpl) ) THEN
          IF ( calc_strair .OR. calc_Tsfc ) THEN
             CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
          ENDIF
-      ELSEIF (nsbc == 4) THEN
+      ELSEIF (ksbc == jp_core) THEN
          IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
             CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
@@ -190,5 +209,5 @@
 
       CALL cice2nemo(aice,fr_i, 'T', 1. )
-      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
+      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_cpl) ) THEN
          DO jl=1,ncat
             CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
@@ -218,9 +237,44 @@
          snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
       ENDIF
-      IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart : 
-         &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh belew sea-ice area
-         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
-         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
-         !
+      IF( .NOT. ln_rstart ) THEN
+         IF( nn_ice_embd == 2 ) THEN            ! full embedment (case 2) deplete the initial ssh below sea-ice area
+            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
+            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
+#if defined key_vvl            
+           ! key_vvl necessary? clem: yes for compilation purpose
+            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
+               fse3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+               fse3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+            ENDDO
+            fse3t_a(:,:,:) = fse3t_b(:,:,:)
+            ! Reconstruction of all vertical scale factors at now and before time
+            ! steps
+            ! =============================================================================
+            ! Horizontal scale factor interpolations
+            ! --------------------------------------
+            CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
+            CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
+            CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
+            CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
+            CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+            ! Vertical scale factor interpolations
+            ! ------------------------------------
+            CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
+            CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
+            CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+            CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
+            CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
+            ! t- and w- points depth
+            ! ----------------------
+            fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+            fsdepw_n(:,:,1) = 0.0_wp
+            fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
+            DO jk = 2, jpk
+               fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
+               fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
+               fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
+            END DO
+#endif
+         ENDIF
       ENDIF
  
@@ -232,5 +286,5 @@
 
    
-   SUBROUTINE cice_sbc_in (kt, nsbc)
+   SUBROUTINE cice_sbc_in (kt, ksbc)
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE cice_sbc_in  ***
@@ -238,5 +292,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt   ! ocean time step
-      INTEGER, INTENT(in   ) ::   nsbc ! surface forcing type
+      INTEGER, INTENT(in   ) ::   ksbc ! surface forcing type
 
       INTEGER  ::   ji, jj, jl                   ! dummy loop indices      
@@ -262,5 +316,5 @@
 ! forced and coupled case 
 
-      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
+      IF ( (ksbc == jp_flx).OR.(ksbc == jp_cpl) ) THEN
 
          ztmpn(:,:,:)=0.0
@@ -287,5 +341,5 @@
 
 ! Surface downward latent heat flux (CI_5)
-         IF (nsbc == 2) THEN
+         IF (ksbc == jp_flx) THEN
             DO jl=1,ncat
                ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
@@ -316,5 +370,5 @@
 ! GBM conductive flux through ice (CI_6)
 !  Convert to GBM
-            IF (nsbc == 2) THEN
+            IF (ksbc == jp_flx) THEN
                ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
             ELSE
@@ -325,5 +379,5 @@
 ! GBM surface heat flux (CI_7)
 !  Convert to GBM
-            IF (nsbc == 2) THEN
+            IF (ksbc == jp_flx) THEN
                ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 
             ELSE
@@ -333,5 +387,5 @@
          ENDDO
 
-      ELSE IF (nsbc == 4) THEN
+      ELSE IF (ksbc == jp_core) THEN
 
 ! Pass CORE forcing fields to CICE (which will calculate heat fluxes etc itself)
@@ -375,9 +429,11 @@
 
 ! Snowfall
-! Ensure fsnow is positive (as in CICE routine prepare_forcing)  
+! Ensure fsnow is positive (as in CICE routine prepare_forcing)
+      IF( iom_use('snowpre') )   CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit  
       ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)  
       CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
 
 ! Rainfall
+      IF( iom_use('precip') )   CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
       ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
       CALL nemo2cice(ztmp,frain,'T', 1. ) 
@@ -458,5 +514,5 @@
 
 
-   SUBROUTINE cice_sbc_out (kt,nsbc)
+   SUBROUTINE cice_sbc_out (kt,ksbc)
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE cice_sbc_out  ***
@@ -464,5 +520,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT( in  ) ::   kt   ! ocean time step
-      INTEGER, INTENT( in  ) ::   nsbc ! surface forcing type
+      INTEGER, INTENT( in  ) ::   ksbc ! surface forcing type
       
       INTEGER  ::   ji, jj, jl                 ! dummy loop indices
@@ -510,5 +566,5 @@
 ! Freshwater fluxes 
 
-      IF (nsbc == 2) THEN
+      IF (ksbc == jp_flx) THEN
 ! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip)
 ! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below
@@ -516,7 +572,7 @@
 ! Better to use evap and tprecip? (but for now don't read in evap in this case)
          emp(:,:)  = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
-      ELSE IF (nsbc == 4) THEN
+      ELSE IF (ksbc == jp_core) THEN
          emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)        
-      ELSE IF (nsbc ==5) THEN
+      ELSE IF (ksbc == jp_cpl) THEN
 ! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above) 
 ! This is currently as required with the coupling fields from the UM atmosphere
@@ -524,6 +580,11 @@
       ENDIF
 
+#if defined key_cice4
       CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
       CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
+#else
+      CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
+      CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
+#endif
 
 ! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
@@ -535,5 +596,6 @@
       sfx(:,:)=ztmp2(:,:)*1000.0
       emp(:,:)=emp(:,:)-ztmp1(:,:)
- 
+      fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
+      
       CALL lbc_lnk( emp , 'T', 1. )
       CALL lbc_lnk( sfx , 'T', 1. )
@@ -543,10 +605,10 @@
 ! Scale qsr and qns according to ice fraction (bulk formulae only)
 
-      IF (nsbc == 4) THEN
+      IF (ksbc == jp_core) THEN
          qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
          qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
       ENDIF
 ! Take into account snow melting except for fully coupled when already in qns_tot
-      IF (nsbc == 5) THEN
+      IF (ksbc == jp_cpl) THEN
          qsr(:,:)= qsr_tot(:,:)
          qns(:,:)= qns_tot(:,:)
@@ -557,5 +619,9 @@
 ! Now add in ice / snow related terms
 ! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
+#if defined key_cice4
       CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
+#else
+      CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
+#endif
       qsr(:,:)=qsr(:,:)+ztmp1(:,:)
       CALL lbc_lnk( qsr , 'T', 1. )
@@ -567,5 +633,9 @@
       ENDDO
 
+#if defined key_cice4
       CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
+#else
+      CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
+#endif
       qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
 
@@ -575,5 +645,5 @@
 
       CALL cice2nemo(aice,fr_i,'T', 1. )
-      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
+      IF ( (ksbc == jp_flx).OR.(ksbc == jp_cpl) ) THEN
          DO jl=1,ncat
             CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
@@ -611,5 +681,4 @@
 
 
-#if defined key_oasis3 || defined key_oasis4
    SUBROUTINE cice_sbc_hadgam( kt )
       !!---------------------------------------------------------------------
@@ -653,10 +722,4 @@
    END SUBROUTINE cice_sbc_hadgam
 
-#else
-   SUBROUTINE cice_sbc_hadgam( kt )    ! Dummy routine
-      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
-      WRITE(*,*) 'cice_sbc_hadgam: You should not have seen this print! error?'
-   END SUBROUTINE cice_sbc_hadgam
-#endif
 
    SUBROUTINE cice_sbc_final
@@ -713,4 +776,25 @@
       IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
          !                                      ! ====================== !
+         ! namsbc_cice is not yet in the reference namelist
+         ! set file information (default values)
+         cn_dir = './'       ! directory in which the model is executed
+
+         ! (NB: frequency positive => hours, negative => months)
+         !            !    file          ! frequency !  variable    ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! landmask
+         !            !    name          !  (hours)  !   name       !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      ! file
+         sn_snow = FLD_N( 'snowfall_1m'  ,    -1.    ,  'snowfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    ) 
+         sn_rain = FLD_N( 'rainfall_1m'  ,    -1.    ,  'rainfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    ) 
+         sn_sblm = FLD_N( 'sublim_1m'    ,    -1.    ,  'sublim'    ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_top1 = FLD_N( 'topmeltn1_1m' ,    -1.    ,  'topmeltn1' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_top2 = FLD_N( 'topmeltn2_1m' ,    -1.    ,  'topmeltn2' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_top3 = FLD_N( 'topmeltn3_1m' ,    -1.    ,  'topmeltn3' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_top4 = FLD_N( 'topmeltn4_1m' ,    -1.    ,  'topmeltn4' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_top5 = FLD_N( 'topmeltn5_1m' ,    -1.    ,  'topmeltn5' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_bot1 = FLD_N( 'botmeltn1_1m' ,    -1.    ,  'botmeltn1' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_bot2 = FLD_N( 'botmeltn2_1m' ,    -1.    ,  'botmeltn2' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_bot3 = FLD_N( 'botmeltn3_1m' ,    -1.    ,  'botmeltn3' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_bot4 = FLD_N( 'botmeltn4_1m' ,    -1.    ,  'botmeltn4' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_bot5 = FLD_N( 'botmeltn5_1m' ,    -1.    ,  'botmeltn5' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+
          REWIND( numnam_ref )              ! Namelist namsbc_cice in reference namelist : 
          READ  ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
@@ -1001,9 +1085,9 @@
 CONTAINS
 
-   SUBROUTINE sbc_ice_cice ( kt, nsbc )     ! Dummy routine
+   SUBROUTINE sbc_ice_cice ( kt, ksbc )     ! Dummy routine
       WRITE(*,*) 'sbc_ice_cice: You should not have seen this print! error?', kt
    END SUBROUTINE sbc_ice_cice
 
-   SUBROUTINE cice_sbc_init (nsbc)    ! Dummy routine
+   SUBROUTINE cice_sbc_init (ksbc)    ! Dummy routine
       WRITE(*,*) 'cice_sbc_init: You should not have seen this print! error?'
    END SUBROUTINE cice_sbc_init
