Index: LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 4523)
+++ LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 5022)
@@ -22,5 +22,10 @@
        radsol, snow, agesno, & 
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa)
+       tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa &
+#ifdef ISO
+       ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
+       xtsnow,xtevap,h1 &  
+#endif            
+       )
 !
 ! This subroutine treats the "open ocean", all grid points that are not entierly covered
@@ -36,4 +41,13 @@
     USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
     use config_ocean_skin_m, only: activate_ocean_skin
+#ifdef ISO
+    USE infotrac_phy, ONLY: ntiso,niso
+    USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall    
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
+    USE isotopes_verif_mod
+#endif
+#endif
 
     INCLUDE "YOMCST.h"
@@ -57,4 +71,10 @@
     real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
 
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(IN)  :: xtprecip_rain, xtprecip_snow
+    REAL, DIMENSION(ntiso,klon), INTENT(IN)  :: xtspechum
+    REAL, DIMENSION(klon),       INTENT(IN)  :: rlat
+#endif
+
 ! In/Output arguments
 !****************************************************************************************
@@ -62,5 +82,9 @@
     REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
     REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno !? put to 0 in ocean
-  
+#ifdef ISO     
+    REAL, DIMENSION(niso,klon), INTENT(IN)   :: xtsnow
+    REAL, DIMENSION(niso,klon), INTENT(INOUT):: Roce
+#endif
+
 ! Output arguments
 !****************************************************************************************
@@ -72,4 +96,9 @@
     REAL, intent(out):: sens_prec_liq(:) ! (knon)
 
+#ifdef ISO     
+    REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux
+    REAL, DIMENSION(klon),       INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation
+#endif
+
 ! Local variables
 !****************************************************************************************
@@ -80,6 +109,10 @@
     REAL, DIMENSION(klon)       :: u1_lay, v1_lay
     LOGICAL                     :: check=.FALSE.
-    REAL sens_prec_sol(knon) 
-    REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol    
+    REAL, DIMENSION(knon)       :: sens_prec_sol
+    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    
+#ifdef ISO   
+    REAL, PARAMETER :: t_coup = 273.15      
+#endif
+
 
 !****************************************************************************************
@@ -87,5 +120,20 @@
 !****************************************************************************************
     IF (check) WRITE(*,*)' Entering ocean_forced_noice'
-    
+
+#ifdef ISO
+#ifdef ISOVERIF
+    DO i = 1, knon
+      IF (iso_eau > 0) THEN         
+        CALL iso_verif_egalite_choix(xtspechum(iso_eau,i), &
+     &                  spechum(i),'ocean_forced_mod 111', &
+     &                  errmax,errmaxrel)     
+        CALL iso_verif_egalite_choix(snow(i), &
+     &                  xtsnow(iso_eau,i),'ocean_forced_mod 117', &
+     &                  errmax,errmaxrel)
+      ENDIF !IF (iso_eau > 0) THEN
+    ENDDO !DO i=1,knon
+#endif      
+#endif 
+
 !****************************************************************************************
 ! 1)    
@@ -103,5 +151,9 @@
 
     else ! GCM
-      CALL limit_read_sst(knon,knindex,tsurf_lim)
+      CALL limit_read_sst(knon,knindex,tsurf_lim &
+#ifdef ISO
+     &     ,Roce,rlat &
+#endif     
+     &     )
     endif ! knon
 !sb--
@@ -161,4 +213,27 @@
          flux_u1, flux_v1)  
 
+#ifdef ISO     
+    CALL calcul_iso_surf_oce_vectall(klon, knon,t_coup, &
+     &    ps,tsurf_new,spechum,u1_lay, v1_lay, xtspechum, &
+     &    evap, Roce,xtevap,h1 &
+#ifdef ISOTRAC
+     &    ,knindex &
+#endif
+     &    )
+#endif         
+
+#ifdef ISO
+#ifdef ISOVERIF
+!          write(*,*) 'ocean_forced_mod 176: sortie de ocean_forced_noice'
+    IF (iso_eau > 0) THEN
+      DO i = 1, knon               
+        CALL iso_verif_egalite_choix(snow(i), &
+     &          xtsnow(iso_eau,i),'ocean_forced_mod 180', &
+     &          errmax,errmaxrel)
+      ENDDO ! DO j=1,knon
+    ENDIF !IF (iso_eau > 0) THEN
+#endif
+#endif   
+
   END SUBROUTINE ocean_forced_noice
 !
@@ -173,5 +248,10 @@
        radsol, snow, qsol, agesno, tsoil, &
        qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l, rhoa)
+       tsurf_new, dflux_s, dflux_l, rhoa &
+#ifdef ISO
+       ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &
+       xtsnow, xtsol,xtevap,Rland_ice &  
+#endif            
+       )
 !
 ! This subroutine treats the ocean where there is ice. 
@@ -187,4 +267,13 @@
     USE indice_sol_mod
     USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
+#ifdef ISO
+    USE infotrac_phy, ONLY: niso, ntiso
+    USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall 
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
+    USE isotopes_verif_mod
+#endif
+#endif
 
 !   INCLUDE "indicesol.h"
@@ -209,4 +298,10 @@
     REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
     real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
+    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum
+    REAL, DIMENSION(niso,klon),  INTENT(IN) :: Roce 
+    REAL, DIMENSION(niso,klon),  INTENT(IN) :: Rland_ice
+#endif
 
 ! In/Output arguments
@@ -216,4 +311,8 @@
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+#ifdef ISO     
+    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: xtsnow
+    REAL, DIMENSION(niso,klon), INTENT(IN)        :: xtsol
+#endif
 
 ! Output arguments
@@ -226,4 +325,7 @@
     REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
+#ifdef ISO     
+    REAL, DIMENSION(ntiso,klon), INTENT(OUT)      :: xtevap
+#endif      
 
 ! Local variables
@@ -238,7 +340,21 @@
     REAL, DIMENSION(klon)       :: u0, v0
     REAL, DIMENSION(klon)       :: u1_lay, v1_lay
-    REAL sens_prec_liq(knon), sens_prec_sol (knon)   
+    REAL, DIMENSION(knon)       :: sens_prec_liq, sens_prec_sol
     REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    
 
+#ifdef ISO
+    REAL, PARAMETER :: t_coup = 273.15
+    REAL, DIMENSION(klon) :: fq_fonte_diag
+    REAL, DIMENSION(klon) :: fqfonte_diag
+    REAL, DIMENSION(klon) :: snow_evap_diag 
+    REAL, DIMENSION(klon) :: fqcalving_diag 
+    REAL, DIMENSION(klon) :: run_off_lic_diag
+    REAL :: coeff_rel_diag
+    REAL :: max_eau_sol_diag  
+    REAL, DIMENSION(klon) :: runoff_diag    
+    INTEGER IXT
+    REAL, DIMENSION(niso,klon) :: xtsnow_prec, xtsol_prec
+    REAL, DIMENSION(klon) :: snow_prec, qsol_prec  
+#endif
 
 !****************************************************************************************
@@ -307,7 +423,68 @@
 !
 !****************************************************************************************
+#ifdef ISO
+   ! verif
+#ifdef ISOVERIF
+    DO i = 1, knon
+      IF (iso_eau > 0) THEN
+        IF (snow(i) > ridicule) THEN
+          CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
+   &              'interfsurf 964',errmax,errmaxrel)
+        ENDIF !IF ((snow(i) > ridicule)) THEN
+      ENDIF !IF (iso_eau > 0) THEN      
+    ENDDO !DO i=1,knon  
+#endif
+   ! end verif
+
+    DO i = 1, knon
+      snow_prec(i) = snow(i)
+      DO ixt = 1, niso
+      xtsnow_prec(ixt,i) = xtsnow(ixt,i)
+      ENDDO !DO ixt=1,niso
+      ! initialisation:
+      fq_fonte_diag(i) = 0.0
+      fqfonte_diag(i)  = 0.0
+      snow_evap_diag(i)= 0.0
+    ENDDO !DO i=1,knon
+#endif
+
+
     CALL fonte_neige( knon, is_sic, knindex, dtime, &
          tsurf_tmp, precip_rain, precip_snow, &
-         snow, qsol, tsurf_new, evap)
+         snow, qsol, tsurf_new, evap &
+#ifdef ISO    
+     &  ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
+     &  ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag   &
+#endif
+     &   )
+
+
+#ifdef ISO
+! isotopes: tout est externalisé
+!#ifdef ISOVERIF
+!        write(*,*) 'ocean_forced_mod 377: call calcul_iso_surf_sic_vectall'
+!        write(*,*) 'klon,knon=',klon,knon
+!#endif
+    CALL calcul_iso_surf_sic_vectall(klon,knon, &
+     &          evap,snow_evap_diag,Tsurf_new,Roce,snow, &
+     &          fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
+     &          precip_snow,xtprecip_snow,xtprecip_rain, snow_prec,xtsnow_prec, &
+     &          xtspechum,spechum,ps, &
+     &          xtevap,xtsnow,fqcalving_diag, &
+     &          knindex,is_sic,run_off_lic_diag,coeff_rel_diag,Rland_ice &
+     &   )
+#ifdef ISOVERIF
+        !write(*,*) 'ocean_forced_mod 391: sortie calcul_iso_surf_sic_vectall'
+    IF (iso_eau > 0) THEN
+      DO i = 1, knon  
+        CALL iso_verif_egalite_choix(snow(i), &
+     &           xtsnow(iso_eau,i),'ocean_forced_mod 396', &
+     &           errmax,errmaxrel)
+      ENDDO ! DO j=1,knon
+    ENDIF !IF (iso_eau > 0) then
+#endif 
+!#ifdef ISOVERIF
+#endif   
+!#ifdef ISO
     
 ! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
