Index: trunk/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limthd_2.F90	(revision 1156)
+++ trunk/NEMO/LIM_SRC_2/limthd_2.F90	(revision 1218)
@@ -7,4 +7,5 @@
    !!            2.0  !  02-07 (C. Ethe, G. Madec) F90
    !!            2.0  !  03-08 (C. Ethe)  add lim_thd_init
+   !!             -   !  08-2008  (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface
    !!---------------------------------------------------------------------
 #if defined key_lim2
@@ -15,5 +16,4 @@
    !!   lim_thd_init_2 : initialisation of sea-ice thermodynamic
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE phycst          ! physical constants
    USE dom_oce         ! ocean space and time domain variables
@@ -31,4 +31,5 @@
    USE limtab_2
    USE prtctl          ! Print control
+   USE cpl_oasis3, ONLY : lk_cpl
       
    IMPLICIT NONE
@@ -37,9 +38,9 @@
    PUBLIC   lim_thd_2  ! called by lim_step
 
-   REAL(wp)  ::   epsi20 = 1.e-20   ,  &  ! constant values
-      &           epsi16 = 1.e-16   ,  &
-      &           epsi04 = 1.e-04   ,  &
-      &           zzero  = 0.e0     ,  &
-      &           zone   = 1.e0
+   REAL(wp) ::   epsi20 = 1.e-20   ! constant values
+   REAL(wp) ::   epsi16 = 1.e-16   !
+   REAL(wp) ::   epsi04 = 1.e-04   !
+   REAL(wp) ::   rzero  = 0.e0     !
+   REAL(wp) ::   rone   = 1.e0     !
 
    !! * Substitutions
@@ -47,5 +48,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!-------- -------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM 2.0,  UCL-LOCEAN-IPSL (2008) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -74,21 +75,20 @@
       INTEGER, INTENT(in) ::   kt     ! number of iteration
       !!
-      INTEGER  ::   ji, jj,    &   ! dummy loop indices
-         nbpb  ,               &   ! nb of icy pts for thermo. cal.
-         nbpac                     ! nb of pts for lateral accretion 
+      INTEGER  ::   ji, jj               ! dummy loop indices
+      INTEGER  ::   nbpb                 ! nb of icy pts for thermo. cal.
+      INTEGER  ::   nbpac                ! nb of pts for lateral accretion 
       CHARACTER (len=22) :: charout
-      REAL(wp) ::  &
-         zfric_umin = 5e-03 ,  &   ! lower bound for the friction velocity
-         zfric_umax = 2e-02        ! upper bound for the friction velocity
-      REAL(wp) ::   &
-         zinda              ,  &   ! switch for test. the val. of concen.
-         zindb, zindg       ,  &   ! switches for test. the val of arg
-         za , zh, zthsnice  ,  &
-         zfric_u            ,  &   ! friction velocity 
-         zfnsol             ,  &   ! total non solar heat
-         zfontn             ,  &   ! heat flux from snow thickness
-         zfntlat, zpareff          ! test. the val. of lead heat budget
-      REAL(wp), DIMENSION(jpi,jpj) ::   zhicifp,   &  ! ice thickness for outputs
-         &                              zqlbsbq       ! link with lead energy budget qldif
+      REAL(wp) ::   zfric_umin = 5e-03   ! lower bound for the friction velocity
+      REAL(wp) ::   zfric_umax = 2e-02   ! upper bound for the friction velocity
+      REAL(wp) ::   zinda                ! switch for test. the val. of concen.
+      REAL(wp) ::   zindb, zindg         ! switches for test. the val of arg
+      REAL(wp) ::   za , zh, zthsnice    !
+      REAL(wp) ::   zfric_u              ! friction velocity 
+      REAL(wp) ::   zfnsol               ! total non solar heat
+      REAL(wp) ::   zfontn               ! heat flux from snow thickness
+      REAL(wp) ::   zfntlat, zpareff     ! test. the val. of lead heat budget
+      REAL(wp) ::   zfi                  ! temporary scalar
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zhicifp   ! ice thickness for outputs
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zqlbsbq   ! link with lead energy budget qldif
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmsk      ! working array
       !!-------------------------------------------------------------------
@@ -100,5 +100,5 @@
       !-------------------------------------------!
       
-!i est-ce utile?  oui au moins en partie
+!!gm needed?  yes at least for some of these arrays 
       rdvosif(:,:) = 0.e0   ! variation of ice volume at surface
       rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom
@@ -114,11 +114,15 @@
       zmsk (:,:,:) = 0.e0
 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            hsnif(ji,jj)  = hsnif(ji,jj) *  MAX( zzero, SIGN( zone , hsnif(ji,jj) - epsi04 ) )
-         END DO
-      END DO
-
-      IF(ln_ctl)   CALL prt_ctl(tab2d_1=hsnif     , clinfo1=' lim_thd: hsnif   : ')
+      ! set to zero snow thickness smaller than epsi04
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            hsnif(ji,jj)  = hsnif(ji,jj) *  MAX( rzero, SIGN( rone , hsnif(ji,jj) - epsi04 ) )
+         END DO
+      END DO
+!!gm better coded (do not use SIGN...)
+!     WHERE( hsnif(:,:) < epsi04 )   hsnif(:,:) = 0.e0
+!!gm
+
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=hsnif, clinfo1=' lim_thd: hsnif   : ' )
       
       !-----------------------------------!
@@ -129,20 +133,20 @@
          DO ji = 1, jpi
             !  snow is transformed into ice if the original ice cover disappears.
-            zindg         = tms(ji,jj) *  MAX( zzero , SIGN( zone , -hicif(ji,jj) ) )
+            zindg         = tms(ji,jj) *  MAX( rzero , SIGN( rone , -hicif(ji,jj) ) )
             hicif(ji,jj)  = hicif(ji,jj) + zindg * rhosn * hsnif(ji,jj) / rau0
-            hsnif(ji,jj)  = ( zone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn
+            hsnif(ji,jj)  = ( rone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn
             dmgwi(ji,jj)  = zindg * (1.0 - frld(ji,jj)) * rhoic * hicif(ji,jj)   ! snow/ice mass
             
             !  the lead fraction, frld, must be little than or equal to amax (ice ridging).
             zthsnice      = hsnif(ji,jj) + hicif(ji,jj)
-            zindb         = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) ) 
-            za            = zindb * MIN( zone, ( 1.0 - frld(ji,jj) ) * uscomi )
+            zindb         = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) ) 
+            za            = zindb * MIN( rone, ( 1.0 - frld(ji,jj) ) * uscomi )
             hsnif (ji,jj) = hsnif(ji,jj)  * za
             hicif (ji,jj) = hicif(ji,jj)  * za
             qstoif(ji,jj) = qstoif(ji,jj) * za
-            frld  (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za , epsi20 )
+            frld  (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za, epsi20 )
             
             !  the in situ ice thickness, hicif, must be equal to or greater than hiclim.
-            zh            = MAX( zone , zindb * hiclim  / MAX( hicif(ji,jj) , epsi20 ) )
+            zh            = MAX( rone , zindb * hiclim  / MAX( hicif(ji,jj), epsi20 ) )
             hsnif (ji,jj) = hsnif(ji,jj)  * zh
             hicif (ji,jj) = hicif(ji,jj)  * zh
@@ -153,9 +157,9 @@
 
       IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=hicif  , clinfo1=' lim_thd: hicif   : ')
-         CALL prt_ctl(tab2d_1=hsnif  , clinfo1=' lim_thd: hsnif   : ')
-         CALL prt_ctl(tab2d_1=dmgwi  , clinfo1=' lim_thd: dmgwi   : ')
-         CALL prt_ctl(tab2d_1=qstoif , clinfo1=' lim_thd: qstoif  : ')
-         CALL prt_ctl(tab2d_1=frld   , clinfo1=' lim_thd: frld    : ')
+         CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif   : ' )
+         CALL prt_ctl( tab2d_1=hsnif , clinfo1=' lim_thd: hsnif   : ' )
+         CALL prt_ctl( tab2d_1=dmgwi , clinfo1=' lim_thd: dmgwi   : ' )
+         CALL prt_ctl( tab2d_1=qstoif, clinfo1=' lim_thd: qstoif  : ' )
+         CALL prt_ctl( tab2d_1=frld  , clinfo1=' lim_thd: frld    : ' )
       ENDIF
 
@@ -175,7 +179,7 @@
          DO ji = 1, jpi
             zthsnice       = hsnif(ji,jj) + hicif(ji,jj)
-            zindb          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) ) 
+            zindb          = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) ) 
             pfrld(ji,jj)   = frld(ji,jj)
-            zinda          = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )
+            zinda          = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
             
             !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
@@ -189,12 +193,24 @@
                         
             !  partial computation of the lead energy budget (qldif)
-            zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting
+#if defined key_coupled 
+            zfi = 1.0 - pfrld(ji,jj)
+            qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                             &
+               &    * (   ( qsr_tot(ji,jj) - qsr_ice(ji,jj) * zfi ) * ( 1.0 - thcm(ji,jj) )   &
+               &        + ( qns_tot(ji,jj) - qns_ice(ji,jj) * zfi )                           &
+               &        + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) )   )
+#else
+            zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting solid precipitation
             zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean
             qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   &
                &                               + zfnsol + fdtcn(ji,jj) - zfontn     &
                &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   &
-               &                               * frld(ji,jj) * rdt_ice    
+               &                        * frld(ji,jj) * rdt_ice    
+!!$            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj) 
+!!$               &           * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )      &
+!!$               &             + qns(ji,jj)  + fdtcn(ji,jj) - zfontn     &
+!!$               &             + ( 1.0 - zindb ) * fsbbq(ji,jj)      )   &
+#endif
             !  parlat : percentage of energy used for lateral ablation (0.0) 
-            zfntlat        = 1.0 - MAX( zzero , SIGN( zone ,  - qldif(ji,jj) ) )
+            zfntlat        = 1.0 - MAX( rzero , SIGN( rone ,  - qldif(ji,jj) ) )
             zpareff        = 1.0 + ( parlat - 1.0 ) * zinda * zfntlat
             zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( (1.0 - frld(ji,jj)) * rdt_ice , epsi16 )
@@ -243,6 +259,6 @@
       !------------------------------------------------------------------------------------ 
 
-      IF ( nbpb > 0) THEN
-         
+      IF( nbpb > 0 ) THEN
+         !    
          !  put the variable in a 1-D array for thermodynamics process
          CALL tab_2d_1d_2( nbpb, frld_1d    (1:nbpb)     , frld       , jpi, jpj, npb(1:nbpb) )
@@ -257,8 +273,8 @@
          CALL tab_2d_1d_2( nbpb, fr2_i0_1d  (1:nbpb)     , fr2_i0     , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, qns_ice_1d (1:nbpb)     , qns_ice    , jpi, jpj, npb(1:nbpb) )
-#if ! defined key_coupled
-         CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb)     , qla_ice    , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d_2( nbpb, dqla_ice_1d(1:nbpb)     , dqla_ice   , jpi, jpj, npb(1:nbpb) )
-#endif
+         IF( .NOT. lk_cpl ) THEN 
+            CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb)     , qla_ice    , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d_2( nbpb, dqla_ice_1d(1:nbpb)     , dqla_ice   , jpi, jpj, npb(1:nbpb) )
+         ENDIF
          CALL tab_2d_1d_2( nbpb, dqns_ice_1d(1:nbpb)     , dqns_ice   , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, tfu_1d     (1:nbpb)     , tfu        , jpi, jpj, npb(1:nbpb) )
@@ -271,7 +287,7 @@
          CALL tab_2d_1d_2( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) )
- 
+         !
          CALL lim_thd_zdf_2( 1, nbpb )       !  compute ice growth
-         
+         !
          !  back to the geographic grid.
          CALL tab_1d_2d_2( nbpb, frld       , npb, frld_1d   (1:nbpb)     , jpi, jpj )
@@ -295,21 +311,20 @@
          CALL tab_1d_2d_2( nbpb, fdvolif    , npb, dvlbq_1d  (1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, rdvonif    , npb, dvnbq_1d  (1:nbpb)     , jpi, jpj ) 
-
- 
-      ENDIF
-
-      
-      !      Up-date sea ice thickness.
-      !---------------------------------
+         !
+      ENDIF
+
+      
+      ! Up-date sea ice thickness
+      !--------------------------
       DO jj = 1, jpj
          DO ji = 1, jpi
             phicif(ji,jj) = hicif(ji,jj)  
-            hicif(ji,jj)  = hicif(ji,jj) *  ( zone -  MAX( zzero, SIGN( zone, - ( 1.0 - frld(ji,jj) ) ) ) )
-         END DO
-      END DO
-
-      
-      !      Tricky trick : add 2 to frld in the Southern Hemisphere.
-      !----------------------------------------------------------
+            hicif(ji,jj)  = hicif(ji,jj) *  ( rone -  MAX( rzero, SIGN( rone, - ( 1.0 - frld(ji,jj) ) ) ) )
+         END DO
+      END DO
+
+      
+      ! Tricky trick : add 2 to frld in the Southern Hemisphere
+      !--------------------------------------------------------
       IF( fcor(1,1) < 0.e0 ) THEN
          DO jj = 1, njeqm1
@@ -321,6 +336,6 @@
       
       
-      !     Select points for lateral accretion (this occurs when heat exchange
-      !     between ice and ocean is negative; ocean losing heat) 
+      ! Select points for lateral accretion (this occurs when heat exchange
+      ! between ice and ocean is negative; ocean losing heat) 
       !-----------------------------------------------------------------
       nbpac = 0
@@ -341,11 +356,9 @@
       ENDIF
 
-      
-      !
-      !     If ocean gains heat do nothing ; otherwise, one performs lateral accretion
+
+      ! If ocean gains heat do nothing ; otherwise, one performs lateral accretion
       !--------------------------------------------------------------------------------
-
       IF( nbpac > 0 ) THEN
-         
+         !
          !...Put the variable in a 1-D array for lateral accretion
          CALL tab_2d_1d_2( nbpac, frld_1d   (1:nbpac)     , frld       , jpi, jpj, npac(1:nbpac) )
@@ -361,8 +374,7 @@
          CALL tab_2d_1d_2( nbpac, dvlbq_1d  (1:nbpac)     , fdvolif    , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d_2( nbpac, tfu_1d    (1:nbpac)     , tfu        , jpi, jpj, npac(1:nbpac) )
-        
-         !  call lateral accretion routine.
-         CALL lim_thd_lac_2( 1 , nbpac )
-         
+         !
+         CALL lim_thd_lac_2( 1 , nbpac )         ! lateral accretion routine.
+         !
          !   back to the geographic grid
          CALL tab_1d_2d_2( nbpac, frld       , npac(1:nbpac), frld_1d   (1:nbpac)     , jpi, jpj )
@@ -375,12 +387,11 @@
          CALL tab_1d_2d_2( nbpac, rdmicif    , npac(1:nbpac), rdmicif_1d(1:nbpac)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpac, fdvolif    , npac(1:nbpac), dvlbq_1d  (1:nbpac)     , jpi, jpj )
-        
+         !
       ENDIF
        
        
-      !      Recover frld values between 0 and 1 in the Southern Hemisphere (tricky trick)
-      !      Update daily thermodynamic ice production.    
+      ! Recover frld values between 0 and 1 in the Southern Hemisphere (tricky trick)
+      ! Update daily thermodynamic ice production.    
       !------------------------------------------------------------------------------
-       
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -392,13 +403,13 @@
       IF(ln_ctl) THEN
          CALL prt_ctl_info(' lim_thd  end  ')
-         CALL prt_ctl(tab2d_1=hicif , clinfo1=' lim_thd: hicif   : ', tab2d_2=hsnif , clinfo2=' hsnif  : ')
-         CALL prt_ctl(tab2d_1=frld  , clinfo1=' lim_thd: frld    : ', tab2d_2=hicifp, clinfo2=' hicifp : ')
-         CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif  : ', tab2d_2=pfrld , clinfo2=' pfrld  : ')
-         CALL prt_ctl(tab2d_1=sist  , clinfo1=' lim_thd: sist    : ')
-         CALL prt_ctl(tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1  : ')
-         CALL prt_ctl(tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2  : ')
-         CALL prt_ctl(tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3  : ')
-         CALL prt_ctl(tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn   : ', tab2d_2=qdtcn , clinfo2=' qdtcn  : ')
-         CALL prt_ctl(tab2d_1=qstoif, clinfo1=' lim_thd: qstoif  : ', tab2d_2=fsbbq , clinfo2=' fsbbq  : ')
+         CALL prt_ctl( tab2d_1=hicif      , clinfo1=' lim_thd: hicif   : ', tab2d_2=hsnif , clinfo2=' hsnif  : ' )
+         CALL prt_ctl( tab2d_1=frld       , clinfo1=' lim_thd: frld    : ', tab2d_2=hicifp, clinfo2=' hicifp : ' )
+         CALL prt_ctl( tab2d_1=phicif     , clinfo1=' lim_thd: phicif  : ', tab2d_2=pfrld , clinfo2=' pfrld  : ' )
+         CALL prt_ctl( tab2d_1=sist       , clinfo1=' lim_thd: sist    : ' )
+         CALL prt_ctl( tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1  : ' )
+         CALL prt_ctl( tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2  : ' )
+         CALL prt_ctl( tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3  : ' )
+         CALL prt_ctl( tab2d_1=fdtcn      , clinfo1=' lim_thd: fdtcn   : ', tab2d_2=qdtcn , clinfo2=' qdtcn  : ' )
+         CALL prt_ctl( tab2d_1=qstoif     , clinfo1=' lim_thd: qstoif  : ', tab2d_2=fsbbq , clinfo2=' fsbbq  : ' )
       ENDIF
        !
@@ -422,11 +433,10 @@
          &                hakdif, hnzst  , thth  , parsub, alphs
       !!-------------------------------------------------------------------
-      
-
-      ! Define the initial parameters
-      ! -------------------------
-      REWIND( numnam_ice )
+      !
+      REWIND( numnam_ice )                  ! read namelist
       READ  ( numnam_ice , namicethd )
-      IF(lwp) THEN
+      IF( lk_cpl .AND. parsub /= 0.0 )   CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )
+      !
+      IF(lwp) THEN                          ! control print
          WRITE(numout,*)
          WRITE(numout,*)'lim_thd_init_2: ice parameters for ice thermodynamic computation '
@@ -437,5 +447,5 @@
          WRITE(numout,*)'       minimum ice thickness                                   hiclim       = ', hiclim 
          WRITE(numout,*)'       maximum lead fraction                                   amax         = ', amax 
-         WRITE(numout,*)'       energy stored in brine pocket (=1) or not (=0)	  swiqst       = ', swiqst 
+         WRITE(numout,*)'       energy stored in brine pocket (=1) or not (=0)          swiqst       = ', swiqst 
          WRITE(numout,*)'       numerical carac. of the scheme for diffusion in ice '
          WRITE(numout,*)'       Cranck-Nicholson (=0.5), implicit (=1), explicit (=0)   sbeta        = ', sbeta
@@ -450,15 +460,15 @@
          WRITE(numout,*)'       coefficient for snow density when snow ice formation    alphs        = ', alphs
       ENDIF
-            
+      !          
       uscomi = 1.0 / ( 1.0 - amax )   ! inverse of minimum lead fraction
       rcdsn = hakdif * rcdsn 
       rcdic = hakdif * rcdic
-      
-      IF ( ( hsndif > 100.e0 ) .OR. ( hicdif > 100.e0 ) ) THEN
+      !
+      IF( hsndif > 100.e0 .OR. hicdif > 100.e0 ) THEN
          cnscg = 0.e0
       ELSE
          cnscg = rcpsn / rcpic   ! ratio  rcpsn/rcpic
       ENDIF
- 
+      !
    END SUBROUTINE lim_thd_init_2
 
